@ 


Advanced 
Engineering 
Research 


Theoretical and scientific-practical jourm 
ISSN 2687-1653 @ 


Mechanics 


DOI 10.23947/2687-1653 


VESUMK-QOTISCGIEG 


Advanced Engineering Research 2022. V. 22, no. 3. ISSN 2687-1653 


Theoretical 
and scientific-practical journal 


Advanced 


Eng ineer ing Published since 1999 


Quarterly 
July-September 2022 


Research 


ISSN 2687-1653 


Vol. 22, no. 3 DOI: 10.23947/2687-1653 


Founder and publisher — Federal State Budgetary Educational Institution of Higher Education 
Don State Technical University (DSTU) 
The journal was known as Vestnik of Don State Technical University (until August 2020) 
Included in the list of peer-reviewed scientific editions where the basic research results of doctoral, candidate’s theses 


should be published (State Commission for Academic Degrees and Titles List) in the following research areas: 


01.02.01 — Analytical Mechanics (Engineering Sciences) 

01.02.04 — Deformable Solid Mechanics (Engineering Sciences) 

01.02.04 — Deformable Solid Mechanics (Physicomathematical Sciences) 

01.02.06 — Dynamics, Strength of Machines, Gear, and Equipment (Engineering Sciences) 

05.02.02 — Engineering Science, Drive Systems and Machine Parts (Engineering Sciences) 

05.02.04 — Machine Friction and Wrear (Engineering Sciences) 

05.02.07 — Technology and Equipment of Mechanical and Physicotechnical Processing (Engineering Sciences) 
05.02.08 — Engineering Technology (Engineering Sciences) 

05.02.10 — Welding, Allied Processes and Technologies (Engineering Sciences) 

05.02.11 — Testing Methods and Diagnosis in Machine Building (Engineering Sciences) 

05.13.11 — Software and Mathematical Support of Machines, Complexes and Computer Networks (Engineering Sciences) 
05.13.17 — Foundations of Information Science (Engineering Sciences) 

05.13.18 — Mathematical Simulation, Numerical Methods and Program Systems (Engineering Sciences) 


The journal is indexed and archived in the Russian Science Citation Index (RSCD), in the International Database 
EBSCO (Academic Search Ultimate Magazines and Journals), Dimensions, CyberLeninka, ROAD, GoogleScholar. 
The journal is a member of Directory of Open Access Journals (DOAJ), Association of Science Editors and Publishers 
(ASEP) and CrossRef. 


Certificate of mass media registration D/I Ne ®C 77 — 78854 of 07.08.2020 is issued by the Federal Service for Supervision 
of Communications, Information Technology, and Mass Media 
The issue is prepared by: 


Inna V. Boyko, Marina P. Smirnova (English version) 


Founder’s, Publisher’s and Printery Address: 
Gagarin Sq. 1, Rostov-on-Don, 344003, Russia. Phone: +7 (863) 2-738-372 


E-mail: vestnik@donstu.ru _ http://vestnik-donstu.ru/ 


© Don State Technical University, 2022 


Advanced Engineering Research 2022. V. 22, no. 3. ISSN 2687-1653 


Editorial Board 
Editor-in-Chief — Besarion Ch. Meskhi, Dr.Sci. (Eng.), professor, Don State Technical University (Russian Federation); 


deputy chief editor — Valery P. Dimitrov, Dr.Sci. (Eng.), professor, Don State Technical University (Russian Federation); 
executive editor — Manana G. Komakhidze, Cand.Sci. (Chemistry), Don State Technical University (Russian Federation); 
executive secretary — Nadezhda A. Shevchenko, Don State Technical University (Russian Federation); 


Sergey M. Aizikovich, Dr.Sci. (Phys.-Math.), professor, Don State Technical University (Russian Federation); 

Kamil S. Akhverdiev, Dr.Sci. (Eng.), professor, Rostov State Transport University (Russian Federation); 

Imad R. Antipas, Cand.Sci. (Eng.), Don State Technical University (Russian Federation); 

Hubert Anysz, PhD (Eng.), assistant professor, Warsaw University of Technology (Republic of Poland); 

Ahilan Appathurai, national junior research fellow, Anna University Chennai (India); 

Gultekin Basmaci, professor (Eng.), Mehmet Akif Ersoy University (Turkey); 

Yuri O. Chernyshev, Dr.Sci. (Eng.), professor, Don State Technical University (Russian Federation); 

Evgenii A. Demekhin, Dr.Sci. (Phys.-Math.), professor, Financial University under the RF Government, Krasnodar branch 
(Russian Federation); 

Oleg V. Dvornikov, Dr.Sci. (Eng.), professor, Belarusian State University (Belarus); 

Karen O. Egiazaryan, Dr.Sci. (Eng.), professor, Tampere University of Technology (Finland); 

Victor A. Eremeey, Dr.Sci. (Phys.-Math.), professor, Southern Scientific Center of RAS (Russian Federation); 
Nikolay E. Galushkin, Dr.Sci. (Eng.), professor, Institute of Service and Business, DSTU branch (Russian Federation); 
LaRoux K. Gillespie, Dr.Sci. (Eng.), professor, President-elect of the Society of Manufacturing Engineers (USA); 

Ali M. Hasan, PhD (Computer Engineering), Al Nahrain University (Baghdad, Iraq); 

Huchang Liao, IET Fellow, BCS Fellow, Sichuan University (China); 

Hamid A. Jalab, PhD (Computer Science & IT), University of Malaya (Malaysia); 

Revaz Z. Kavtaradze, Dr.Sci. (Eng.), professor, Raphiel Dvali Institute of Machine Mechanics (Georgia); 

Janusz Witalis Kozubal, Dr.Sci. (Eng.), Wroclaw Polytechnic University (Republic of Poland); 

Ilya I. Kudish, PhD (Phys.-Math.), Kettering University (USA); 

Victor M. Kureychik, Dr.Sci. (Eng.), professor, Southern Federal University (Russian Federation); 

Geny V. Kuznetzov, Dr.Sci. (Phys.-Math.), professor, Tomsk Polytechnic University (Russian Federation); 

Vladimir I. Lysak, Dr.Sci. (Eng.), professor, Volgograd State Technical University (Russian Federation); 

Vladimir I. Marchuk, Dr.Sci. (Eng.), professor, Institute of Service and Business, DSTU branch (Russian Federation); 
Vladimir M. Mladenovic, Dr.Sci. (Eng.), professor, University of Kragujevac (Serbia); 

Murman A. Mukutadze, Cand.Sci. (Eng.), professor, Tambov State Technical University (Russian Federation); 
Andrey V. Nasedkin, Dr.Sci. (Phys.-Math.), professor, Southern Federal University (Russian Federation); 

Tamaz M. Natriashvili, academician, Raphiel Dvali Institute of Machine Mechanics (Georgia); 

Nguyen Dong Ahn, Dr.Sci. (Phys.-Math.), professor, Academy of Sciences and Technologies of Vietnam (Vietnam); 
Nguyen Xuan Chiem, Dr.Sci. (Eng.), Le Quy Don Technical University (Vietnam); 

Sergey G. Parshin, Dr.Sci. (Eng.), associate professor, St. Petersburg Polytechnic University (Russian Federation); 
Konstantin V. Podmaster’ev, Dr.Sci. (Eng.), professor, Orel State University named after I.S. Turgenev (Russian Federation); 
Roman N. Polyakov, Dr.Sci. (Eng.), associate professor, Orel State University named after I.S. Turgenev (Russian Federation); 
Valentin L. Popov, Dr.Sci. (Phys.-Math.), professor, Berlin University of Technology (Germany); 

Nikolay N. Prokopenko, Dr.Sci. (Eng.), professor, Don State Technical University (Russian Federation); 

José Carlos Quadrado, PhD (Electrical Engineering and Computers), Polytechnic Institute of Porto (Portugal); 
Alexander T. Rybak, Dr.Sci. (Eng.), professor, Don State Technical University (Russian Federation); 

Muzaafer H. Saraéevié, full professor, Novi Pazar International University (Serbia); 

Arestak A. Sarukhanyan, Dr.Sci. (Eng.), professor, National University of Architecture and Construction of Armenia 
(Armenia); 

Vladimir N. Sidorov, Dr.Sci. (Eng.), Russian University of Transport (Russian Federation); 

Arkady N. Solovyev, Dr.Sci. (Phys.-Math.), professor, Don State Technical University (Russian Federation); 
Mezhlum A. Sumbatyan, Dr.Sci. (Phys.-Math.), professor, Southern Federal University (Russian Federation); 
Mikhail A. Tamarkin, Dr.Sci. (Eng.), professor, Don State Technical University (Russian Federation); 

Murat Tezer, professor, Near East University (Turkey); 

Bertram Torsten, Dr.Sci. (Eng.), professor, TU Dortmund University (Germany); 

Vyacheslav G. Tsybulin, Dr.Sci. (Phys.-Math.), associate professor, Southern Federal University (Russian Federation); 
Umid M. Turdaliev, Dr.Sci. (Eng.), professor, Andijan Machine-Building Institute (Uzbekistan); 

Ahmet Uyumaz, Dr.Sci. (Eng.), professor, Mehmet Akif Ersoy University (Turkey); 

Valery N. Varavka, Dr.Sci. (Eng.), professor, Don State Technical University (Russian Federation); 

Igor M. Verner, Cand.Sci. (Eng.), docent, Technion (Israel); 

Sergei A. Voronov, Dr.Sci. (Eng.), associate professor, Russian Foundation of Fundamental Research (Russian Federation); 
Batyr M. Yazyev, Dr.Sci. (Phys.-Math.), professor, Don State Technical University (Russian Federation); 

Vilor L. Zakovorotny, Dr.Sci. (Eng.), professor, Don State Technical University (Russian Federation). 


Advanced Engineering Research 2022. V. 22, no. 3. P. 179 ISSN 2687-1653 


CONTENTS 


MECHANICS 


Petukhov D.S., Adamov A.A., Keller I. E. Selection and Identification of a Model of Elasto- 
Viscoplasticity of the Filled Fluorocomposite according to Free and Constrained Compression Tests................ 
Belonozhko D. G., Korolev I. D., Chernyshev Y. O., Ventsov N. N. Methodology for Determining the 
Probability of Damage to Unmanned Aerial Vehicles as a Result of Exposure to Atmospheric Environ- 
MEM tal PACtOrs ssscsnccssssescssescsiishs sveshdecsoassasstedadecsaceisessodshdnstaessssicnacsdestsseanehsnssstesdpsassescodiedssiatagenssuatetessrsedessanabadss ipaasapsoatbivesas 
Saypulaev G. R., Saypulaev M. R., Merkuryev I. V., Adamov B. I., Garcia R. B. Application of an Inertial 
Sensor Unit for Position Estimation and Motion Control of the Lower-Extremity Powered Exoskeleton .......... 
Peshkhoev I. M., Sobol B. V. Critical Loads of Uniformly Compressed Orthotropic Rectangular Plate on 
AN HL AStie Bases. ssescscscsziiess sschecthscabives anssdsseaiseys pontsbeosestessayanseacdodbiess ceecvedsanviseas sgeecosonssevss cueedeovecstéascieesdosesstvesscvoyeabsbetibnsde 
Pakhnutova N. V., Boyangin E. N., Shkoda O. A., Zelepugin S. A. Microhardness and Dynamic Yield 
Strength of Copper Samples upon Impact on a Rigid Wall 


MACHINE BUILDING AND MACHINE SCIENCE 


Lyudmirsky Y. G., Leonov V. P., Assaulenko S. S. Increasing the Durability of Butt-Welded Joints Operat- 
ing under Cyclic Loads in a Biaxial Stress Field .ccccccecccesesssssssssssssssesseensssssssssseeeeesesienesssssssseeeeseesesmmnsssssseeseeeeseeiienensssssseeeeeeee 
Tsybrii I. K., Koval N. S., Issabekov Z. N. Study of Tapered Threads Parameters Affecting the Quality of 
the: Oil Pipelame: COMME CON 5:5: s22ssd5ccseiisssscaheteostecstehsadsstsivszssesconsonsithsidaeseudeadecenstesossnssaliytndisebboossaues eat suncentée; pesrscsdosentsuseses 


Lyudmirsky Y. G., Assaulenko S. S., Kramskoi A. V. Methods and Equipment for Experimental Evaluation 
of the Performance of Shell and Hull Structures 


INFORMATION TECHNOLOGY, COMPUTER SCIENCE, AND MANAGEMENT 


Amera Al-Funjan, Farid Meziane, Rob Aspin. Describing Pulmonary Nodules Using 3D Clustering .............. 
Elshamy M. M. M., Tiraturyan A. N., Uglova E. V., Elgendy M. Z. Evaluation of Pavement Condition 
Deterioration Using Artificial Intelligence Models... eesessssssesceeeceseeeeeeeesesesesearacsnsesucnensssseacacececeseseseeneseseaees 
Babushkina N. E., Lyapin A. A. Solving the Problem of Determining the Mechanical Properties of Road 
Structure Materials Using Neural Network Technologies 


179 


http://vestnik-donstu.ru 


180 


Advanced Engineering Research 2022. V. 22, no. 3. P. 180-192. ISSN 2687-1653 


MECHANICS 


® Check for updates 
BY 


UDC 539.374 Original article 
https://doi.org/10.23947/2687-1653-2022-22-3-180-192 


Selection and Identification of a Model of Elasto-Viscoplasticity of the Filled 
Fluorocomposite according to Free and Constrained Compression Tests 


Dmitriy S. Petukhov , Anatoliy A. Adamov™, Ilya E. Keller 
Institute of Continuous Media Mechanics UrB RAS, 1, Academika Koroleva St., Perm, Russian Federation 


DX] petuhovds @mail.ru 


Abstract 
Introduction. Properties of filled composites based on polytetrafluoroethylene allow them to work at high contact 


pressures, reciprocating nature of shear loads, and in a wide temperature range. Due to this, they are used as antifriction 
layers of bearing parts with ball segment. To simulate the mechanical behavior of such materials under operating 
conditions, adequate constitutive equations of elasto-viscoplasticity and methods of their identification according to the 
data of basic experiments are needed. 

Materials and Methods. The tensor-linear model of elasto-viscoplasticity should be identified according to the data of 
tests on free compression of samples. They were subjected to loading up to a maximum deformation of 10 %, allowed 
to remain, unloaded, and then, a similar loading cycle up to 160 MPa under constrained compression was carried out. 
The experiment with a composite based on polytetrafluoroethylene filled with 40 wt. % fine bronze, was conducted at 


room temperature. Tests for constrained compression were performed for two values of the strain rate, and for free 
compression — for three values of the strain rate in the range of 10°—10~s''. For the description, two models of 


elasto-viscoplasticity were considered, representing modifications of Swain and Kletschkowski’s models and 
corresponding to the connection of a viscoelastic or elastic nonlinear viscous element with a plastic or endochronic 
element. An integral operator with a Kohlrausch kernel was considered as a viscoelastic element. 

Results. The results of the constrained compression tests made it possible to separate the elastic relationship of 
volumetric deformations and average stresses from the constitutive relations. The data of free compression cycles at 
different strain rates were used to determine the material constants of the model. For this purpose, an efficient search 
algorithm based on the simplex method of minimizing the discrepancy was implemented. Both models discovered the 
importance of the plastic component (independent of the deformation rate) for a qualitative description of the stress 
cycling that accompanied the cyclic deformation, as well as their dependence on the strain rate. 

Discussion and Conclusions. Both models of elasto-viscoplasticity described correctly the behavior of the studied 
fluorocomposite under loading conditions close to the operating conditions of the antifriction layers of the bearing parts 
with ball segment. They can be considered as a basis for their further generalization, taking into account the dependence 


on temperature. 


Keywords: elasto-viscoplasticity, constitutive equations, filled fluorocomposite, identification, free and constrained 


compression, trapezoidal loading. 
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Introduction. Properties of filled composites based on polytetrafluoroethylene (PTFE) allow them to work in a 
wide temperature range at high contact pressures and reciprocal sliding on the counterbody. Due to this, they are used as 
antifriction layers of bearing parts with a ball segment'*. PTFE as a base for antifriction polymer composites in a wide 
temperature range shows high stability of properties and exceptionally low coefficients of friction, sliding and rest. The 
introduction of fillers helps to significantly increase the wear resistance, stiffness and yield strength, as well as reduce 
the creep of the composite without marked increase in the friction coefficients. The rheological properties of the 
antifriction material should provide feasibility during the operation of the product in the temperature range from less 


than —50 to +50 °C: 
—at calculated values of the compressive load of at least 60 MPa; 
—at its peak values up to 150 MPa; 
—at a normalized range of cyclic shear deformations. 


The predictability of changes in the properties and thickness of the antifriction layer under its hinge mounting also 


matters. 


During the operation of the materials under consideration, it is required to assess the performance of antifriction 
layers of the support parts with a ball segment, to predict the resource for the entire service life (up to 50 years). To do 
this, you need to construct and identify adequate constitutive equations of elasto-viscoplasticity according to the data of 
basic experiments. Models describing the rheological behavior of filled fluorocomposites are also important for 
calculating the operational properties of shaft seals and rods in high-power diesel engines*. To calculate the behavior of 
products under load during installation (for short periods of time), it is important that the models correctly describe the 
plastic (elastoplastic, viscoplastic) properties of the material. To calculate the behavior of products during operation (for 
long periods of time), it is required to adequately represent the phenomena of relaxation, creep, and cyclic creep 


(ratcheting). 


Rheological models of filled fluorocomposites describing cycles of loading, holding, and unloading in wide ranges 
of time and strain rates are rarely found in the literature. We can specify only work [1] that meets all these requirements. 
It proposes a model corresponding to the parallel connection of endochronic and nonlinear viscous elements. This 


model describes all uniaxial tensile tests* at room temperature. It is about: 


— loading curves up to 0.08 strain with strain rates 10“ —107s"'; 


— relaxation curves for 10 h with constant strain of 0.03—0.08; 


— hysteresis loops in the load cycle up to 0.1 strain, short-term hold and unloading with subsequent hold. 


Mechanics 


' EAD 050009-00-0301. Spherical and cylindrical bearing with special sliding material made of fluoropolymer. EOTA. Product regulation. URL: 
https://www.nlfnorm.cz/en/ehn/6189 

? Stanton JF, Roeder CW, Campbell TI. Appendix C: Friction and Wear of PTFE Sliding Surfaces. NCHRP Report 432: High-Load Multi-Rotational 
Bridge Bearings. Washington, DC: TRB; 1999. 413 p. 

> Hai Sui, Heiko Pohl, Achim Oppermann, et al. Material and Computational Analysis of PTFE Seals. SAE Paper Series. 1995:951055. URL: 
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One-dimensional constitutive equations containing seven material constants were written, as well as their tensor 
generalization for the geometrically linear case. The authors point out that the endochronic plasticity model makes it 
possible to describe a continuous change in plastic deformation during loading and unloading. This is typical for 
polymer materials and is not typical for metals. Note that in [1], the endochronicity parameter was not used, which 
allowed regulating the characteristic range of transition to the plastic state. In addition, the paper’ describes similar 
curves for temperatures up to 120°C, and strain rates in the range of 10*-100 s", relaxation curves with a sharp 


increase in temperature up to 100 °C, and relaxation curves under uniaxial compression. 


A tensor geometrically nonlinear model is presented in [2]. It describes the behavior of filled fluorocomposites 
under various monotonic and cyclic loading histories in various stress states and at arbitrary temperatures. A 
viscoplastic rheological element in a series connection with a viscoelastic one is taken into account. On the whole, this 
model satisfactorily describes the curves corresponding to the experiments [1], but it is much more complicated and 
requires fourteen material constants. In [3], the model [4] is used to describe microindentation cycles of a filled 
fluorocomposite. The rheological block diagram is similar to [2], but it takes into account the relaxation time spectrum 
of the viscoelastic element. In [3], the model [4] is used to describe microindentation cycles of a filled fluorocomposite. 
The rheological block diagram is similar to [2], but it takes into account the relaxation time spectrum of the viscoelastic 
element. 

The model [5] was designed to describe features of the elastoplastic behavior of solid polymers under cyclic uniaxial 
loading. The evolutionary equations for its internal variable, stresses and plastic deformations, had nonlinear cross- 
connections, therefore, the model did not correspond to any rheological structural scheme. The measure of internal time 
characteristic of the endochronic theory of plasticity was used. Twenty-one constants were required to identify the 
model. 

The nonlinear one-dimensional Maxwell-type viscoelastic plasticity model [6-9] contains only two material 
functions. It describes incomplete recovery of deformations in the load-unloading cycle and cyclic creep (ratcheting) 
with asymmetric loading cycles. Nonlinear viscoelasticity models with fractional integro-differentiation operators are 
characterized by flexibility with a small number of material constants and functions [7]. It is necessary to investigate the 
applicability of such nonlinear models to describe the behavior of polymer materials in the loading and unloading cycle 
with holds [10-13]. Papers [14, 15] are devoted to the identification of linear models of viscoelasticity based on the 
results of indentation tests. In [16, 17], the dynamic mechanical analysis results are used for this purpose, but the 
significant nonlinearity of the materials studied in this work does not allow using this technique. Ratcheting of PTFE 


composites is experimentally investigated in [18, 19], but there is no mathematical description of the results. 


Paper [8] presents a method for identifying a model of isotropic elastoplastic behavior of filled composites based on 
PTFE and ultrahigh molecular weight polyethylene based on experimental data on free and constrained compression. 
Due to the peculiarities of the rheology of plastics [9], the methods [10] are applicable to them. This approach makes it 


possible to determine: 
— the elastic volume compressibility function; 
— the hardening function for shear elastoplastic properties within the framework of the theory of plastic flow. 


Thus, the operating conditions of the products in question, the mechanical properties of polymer materials, and the 


presented literature analysis serve as the basis for further research. 


Let us highlight the basic test for free compression of cylindrical samples of 20x20 mm. It consists of load cycles up 
to 0.1 deformation, a short-term hold of fifteen minutes, and unloading with a subsequent hold of fifteen minutes. In this 


way, data can be determined to describe elastic, plastic and viscous components of deformation, relaxation and reverse 


> Material and Computational Analysis of PTFE Seals. 
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creep with residual deformations. A similar definition of hardness is provided by the standard®. Loading is performed at 
three strain rates in the range of 10°-10°s". Besides, the samples are tested for constrained compression to a stress of 
160 MPa at two strain rates from the same range. All tests are repeated for values from the operating temperature range 
of the polymer material. These data should be sufficient to identify a suitable model of elasto-viscoplastic compressible 


material with nonlinear elastic compressibility. 


First of all, we are interested in the possibilities of describing experimental data at room temperature with the 


simplest elasto-viscoplasticity models of the type proposed in [1] and [4], which differ in structural schemes. 


Materials and Methods. A composite based on PTFE PN 90 filled with 40 wt. % of fine bronze was chosen as the 
simulated material. Cylindrical samples with a diameter and height of 20 mm were made by pressing and sintering a 


mixture of powders. The following experiments were conducted to select and then identify the model: 
— for constrained compression with a single nonzero component of the strain tensor; 


— for free compression with a single nonzero component of the stress tensor. 


In the first case, friction along the lateral boundary was excluded by lubrication, in the second — by a thin film of 
PTFE and grease. 


2,000 4,000 6,000 ts 


a) 
o, MPa 


150 


100 


50 


b) 


Fig. 1. Cyclic tight compression test: a — program; b — hysteresis curve 
The experiments were aimed at identifying irreversible volumetric deformations and volume-elastic properties of the 
material. The first test was cyclic constrained compression with increasing amplitude (Fig. 1 a). It set the value of the 
deformation rate €=2.5x10~°s”" and the pause between unloading and loading for 300 s. During the first few cycles, a 


residual volume strain equal to approximately 5x10~° accumulated (Fig. 1). It changed slightly when the sample was held in 


the unloaded state. 


® GOST 4670-2015 (ISO 2039-1:2001). Plastics. Determination of hardness. Ball indentation method. Interstate Council for Standardization, 
Metrology and Certification. Moscow: Standartinform; 2016. 10 p. (In Russ.) 
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Sg and 


Further, compression tests were carried out up to 160 MPa with a significantly slower strain rate € =2.5x10 


subsequent unloading (Fig. 2 a): 
— of undeformed sample (1); 
— of the same sample repeatedly (2); 


—of the sample after the cyclic test shown in Figure 1 (3). 


o, MPa 
150 Test no.1 
— Test no.2 
— Test no.3 
100 
50 
0 
0.01 0.02 0.03 0.04 & 
— Testno.1 
— Test no.2 
— Test no.3 


0.001 0.002 0.003 0.004 € 
b) 


Fig. 2. Slow tight compression tests: a — loading diagrams; b — initial and final sections of diagrams 


Tests of the same sample were carried out in a day. 

As can be seen in Figure 2 b, in the experiment with a new sample, the maximum residual strain of about 3.5x10° is 
observed. On repeated testing, it amounted to 1x10”. In the third test, there is no residual strain. The properties of the 
material under cyclic constrained compression are accompanied by a transient irreversible process and are stabilized in 
several cycles. 

The bulk properties were determined from the third test (Fig. 2 a), in which there are no residual strains. In this 
paper, bulk strains are assumed to be nonlinearly elastic. The data is described by a quadratic dependence: 

o=M(sje; M(e)=a,+a,6; a, =2.4hPa; a, =22.1 hPa. (1) 

Here, 6,6 — axial components of the Cauchy stress tensors and logarithmic strains; M — constrained compression 


modulus, which is related to the bulk compression modulus K_ by the ratio: 


K(e,)(3K(e,) +E 
Mee, )=3 ReIGKE) +B.) (2) 
9K(e,)—-E,, 
where E,, — equilibrium Young's modulus (independent of bulk strain); €, =€;, € — components of the logarithmic 


strain tensor. 
Value E,, =690 MPa was obtained with the slope of the curve o(€) under uniaxial compression in the vicinity of 


o = 1MPa and relaxation 8.65 %. 
Using (1), (2) and 9.38 % correction for the compliance of the machine in the tight compression test, it is possible to 


obtain the desired dependence K(é,) , which differs slightly (maximum by 0.03 %) from linear: 
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0,, =K(e, Jey; K(e,)=b,+b,; by, =2.1hPa; b, =22.2 hPa, (3) 


where 6,, =+6,, 6; — components of the Cauchy stress tensor. 


i> OF 

Experiments to determine the relationship between average stress and volume strain were carried out on samples 
previously tested by constrained compression with irreversible volume strain. The same applies to the free compression 
experiments described below and required to determine the relationship of the deviatory parts of stress and strain 
tensors. As a result, the model reflected the properties of the compacted material. In production, processes associated 
with irreversible bulk strain occur under the crimping of a layer of antifriction polymer material during the assembly of 
the structure of the support part and do not affect the further behavior of the material. 

Uniaxial compression tests were carried out to determine the relationships between the deviatory parts of stress and 


strain tensors. Figure3 shows the results of experiments for three values of strain rates 
&=4.63x10°s', 2.35x10°s', 2.27x10°s''. Further, the data from these experiments will be interpreted in the 


framework of a one-dimensional model. When formulating a three-dimensional model, the elastic volume part can be 
excluded from these relations. 


o, MPa 


Experiments at different speeds 
15 — ¢=4.63x10° 
— €=2.35x10" 
— €=2.27x10" 
10 
5 
0 ZZ | 


0.02 0.04 0.06 008 010 « 
a) 


0.10 Experiments at different speeds 
, — £=4.63x10° 
0.08 - €=2.35x10° 
— €=2.27x10° 
0.06 
0.04 
0.02 


5 10 15 20. 25 1x10°, s 


Experiments at different speeds 
15 - ¢=4.63x10° 
- €=2.35x10° 
— €=2.27x10% 
10 
5 
0 


5 10 15 20 «25 1x10°, s 
c) 


Fig. 3. Uniaxial compression experiment with three different speeds: 
a — loading diagrams; b — deformation histories; c — stress histories 
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When choosing a suitable model, it turned out that the experimental data in Figure 3 were not described by the 
relations of a linear viscoelastic medium: 


o(t)=[R(t-r)e(r)d7 = Ry [ec -{r¢ pena} (4) 


where o(t), €(t) — stress and strain histories used to describe elastomers and polymers [10]. 
There have been attempts to use expressions for relaxation functions R(t) or relaxation rates 7“(t) with a different 


number of parameters: the sum of exponents, the Robotnov fractional-exponential function’, the Koltunov kernel’, the 
kernel from the work of T. L. Smith’, the Kohlrausch kernel [10]: 


—-rR iS ie Thee _ Aexp(-Bt”). 
RO=R, +). R, exp(—a,t); M(t) = At LT (ain) T(t) = : s 
R,-R.,, 7 : 
RD Rots Gaye KO=Re Re Rae tay"). 


Here, R,,, Ry, R,;, @, 8, T — material constants; [ — gamma function, as well as monotonic approximation by 
cubic splines at eight points, whose coordinates served as parameters. 
Models were considered, in which, in addition to viscous ones, plastic irreversible deformations were taken into 


account. Two basic models of this type were investigated. 
Basic model | (Fig. 4a) is a combination of models [4] and [10] and assumes a sequential connection of three 


elements. 


Fig. 4. Block diagrams of basic models: a — model 1; b — model 2 


1. Linear viscoelastic element with Kohlrausch kernel [10]: 


: a 
€,(t) = J T(t—r)6(r)dr, ‘7 (t) = 7 (z em Jos ‘) (6) 
where R, R,,, 7), @ — material constants, the first two have the physical meaning of instantaneous and equilibrium 
modules. 
2. Plastic element with linear hardening: 
. k6(t), C26, Vvo<S-o, 
s10=| 0, -o,<o<o, ° 7) 
where o,,, k — material constants. 
3. Linear-viscous element: 
é,()=1 'o(t), (8) 
where 1 — material constant. 
The final expression for strains: 
e(t) =€,(t) +€,(t)+e,(t). (9) 


7 Rabotnov YuN. Polzuchest' ehlementov konstruktsii. Moscow: Nauka; 1966. 752 p. (In Russ.) 

* Koltunov MA. Singulyarnye funktsii vliyaniya v analize relaksatsionnykh protsessov. Prochnost' i plastichnost'. Moscow: Nauka; 1971. P. 640-645. 
(In Russ). 

° Smith TL. Ehmpiricheskie uravneniya dlya vyazkouprugikh kharakteristik i vychisleniya relaksatsionnykh spektrov. Vyazkouprugaya relaksatsiya v 
polimerakh. Moscow: Mir; 1974. P. 44—56. (In Russ). 
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Model (6)-(9) is solved with respect to strains, therefore, it is convenient to have stress histories for its 
identification o(t) (Fig. 3 c). 

The block diagram of basic model 2 (Fig. 4 b), borrowed from [1], is a parallel connection of two elements. 

1. Endochronic plastic element with nonlinear elastic part: 


é()= SO lac) ; 6,(t) = Asign(e(t) —e,(t))In(1 +Ble(t)—e, (t) ), (10) 
where €, — an internal parameter that has the meaning of irreversible strain of the element, Y,A,B — material 
constants. 

2. A pair of sequentially connected nonlinear-viscous (2) and linear-elastic (3) elements: 
k 
a,(t)= #sin(o, (0) ; (1) 
7] 0 
0,(t) =0,(t) =ce,(t), (12) 
where sign — argument sign; n, 0,, k , c — material constants. 
The resulting expression for stresses: 
o(t) =0,(t)+o,(t). (13) 


Model (10)—(13) is solved with respect to stresses, therefore, it is convenient to have strain histories for its 
identification g(t) (Fig. 3 c). Expression o,(t) in (13) is solved as the solution to a system of algebraic differential 
equations (10), and o,(t) — as a system solution (11), (12) taking into account g(t) =¢,(t)+¢,(t). 

Research Results. Residual minimization between experimental strain histories (Fig. 3 b) and model predictions 
allowed us to find seven material constants of base model 1: 

R =895 MPa;R,, = 205 MPa; t, = 1.66x10's; a = 0.64; 
6, =10.2 MPa; k = 0.0095 MPa; 1 =9.28x10°s. 


The simplex method is used as a search procedure. Comparison of the calculated and experimental data is shown 


(14) 


in Figure 5. 
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Quick experiment 
0.10 =o — Experiment 
— Model 


0.08 e Residual points 


250 500 750 1,000 1,250 1,500 t,s 
c) 
Fig. 5. Comparison of experimental data and base model 1 predictions: 
a— &=4.63x10°s';b— €=2.35x10°s";c— €=2.27x10" s* 


Residual minimization between the experimental stress histories (Fig. 3 b) and the model predictions allowed us to 


detect seven material constants of base model 2: 


Y = 20.3MPa; A = 6.01 MPa; B= 61.8; 


15 
1 =6.81-10° s; 6, = 0.66 MPa; k = 1.08; s = 681 MPa. (15) 


Here, strain histories were considered to be given (Fig. 3 b). Comparison of the calculated and experimental data is 


shown in Figure 6. 
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o, MPa 


Quick experiment 


— Experiment 
— Model 
e Residual points 


250 500 750 1,000 1,250 1,500 tx10°, 


c) 
Fig. 6. Comparison of experimental data and predictions of base model 2: 


a— €=4.63x10°s";b— €=2.35x10°s';c— €=2.27x10" Ss" 
Discussion and Conclusions. The considered one-dimensional models were generalized to a spatial formulation 
that assumed isotropy of properties and immutability of the main axes of the strain tensor under loading. To do this, the 


Cauchy stress tensors and logarithmic deformations were decomposed into spherical and deviatoric parts: 
1 
o=o0,I+58, € mae 


The relationship between the ball parts 6,, and €, is given by expression (3). The relationship between the 


m 


deviatoric parts s and e is given by the relations generalizing (6)—(9) or (10)-(13) depending on the selection of the 
base model. These expressions should be rewritten in terms of the histories of the deviatoric parts s(t) and e(t), 


excluding the ball part from uniaxial compression. 

AS a Spatial generalization of the relations (7) for the plastic element of model 1, the plastic flow law associated with 
von Mises yield criterion, with isotropic linear strain hardening was considered. The connection of deviators will look 
like this: 

e(t) =e,(t) +e, (t) +e, (t); 


e,(t) = | 7(t-1)8(r)dr; M)== (i 1 Jos( (t/t,)*); 


t 


e, =2k ay e,(t) = fé, (od; Gy = fee, 26s; 
e, 0 


-l 
e,(t)=n s(t), 
where colon means the operation of convolution of tensors A: B= A,B, given by the components in an orthonormal 


basis. 
The relations of the endochronic plasticity theory (10) of model 2 are also generalized in the spirit of von Mises 
plasticity, which provides the following entry for deviators: 


S=S,+8,;e=e, +e; 


3 
6, = 56S; &=,/26:6; 5 =dev{ Y'éinn,} é, = Asign (A, )In(1+B 
i=l 


di 


; 


k 


31 fo) S 
2 2 ae =. | 8a. eae & 
é, 5 sh 3 0, =/38, :8); 
n 99 /| S 
=n iy 


where “dev” means a tensor deviator, A,, n; — eigenvalues and eigenvectors of the tensor e—e, =, n, nj. 
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It should be noted that the authors [1] used a model option with an endochronic parameter [14] equal to one. It 
seems appropriate to consider this parameter as a material constant, since its influence on the description of transients 
after changes of loading modes is known. 

In the framework of this work, the experiments were carried out on constrained and free compression of a composite 
based on polytetrafluoroethylene filled with finely dispersed bronze with a mass fraction of 40 %. At the same time, 
loading, unloading and holding were carried out with different strain rates in the range of 10°-10™ s"'. The experiments 
on constrained compression were carried out at pressures corresponding to the peak loads of the antifriction layer in 
hinges with a spherical segment. The inability of the standard linear viscoelastic model to describe test data for several 
strain rates was found. Two families of seven-constant models with plastic and nonlinear viscous structural elements 
that could describe the data of basic experiments were selected and identified. The model with an endochronic element 
showed great accuracy in describing the data at a high strain rate. The model with a plastic element can be improved if 
the nonlinearity in the viscous element and nonlinear strain hardening are taken into account, which will require an 
increase in the number of material constants. The data of these models were generalized to the spatial case for an 


isotropic material and strain histories with a constant orientation of the trihedron of the main axes. 
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Abstract 

Introduction. For the effective use of unmanned aerial vehicles (UAV), it is required to take into account the impact of 
atmospheric environmental factors. Based on the existing level of development of the research and methodological 
apparatus, it is impossible to determine the probability of damage to the UAV under the conditions of complex exposure 
to atmospheric environmental factors and to assess the feasibility of further performance of the flight task. Violation of 
the UAV operation process is caused by the influence of atmospheric precipitation, wind and temperature conditions of 
the environment. The work aims at the development of a methodology for determining the probability of damage to 
UAV as a result of exposure to atmospheric environmental factors, as well as to evaluate the performance of the 
software package implementing the developed algorithm using triangular, trapezoidal, pentagonal, and Gaussian 
membership functions. 

Materials and Methods. A technique is proposed that makes it possible, using the theory of fuzzy logic, to determine 
the probability of damage to the UAV with inaccuracies and uncertainties in the description of atmospheric effects of 
the environment. It takes into account possible atmospheric forcing and enables to determine the probability of damage 
to the UAV under various atmospheric influences. The computational complexity of the algorithm implementing the 
technique depends significantly on the number of qualitative assessments of atmospheric impacts on UAV. 

Results. A method for determining the probability of damage to UAV as a result of exposure to atmospheric 
environmental factors based on fuzzy sets (triangular, trapezoidal, pentagonal, Gaussian) was proposed and tested. The 
use of fuzzy sets for estimating the conditions of the UAV application environment was described. An algorithm was 
developed to determine the probability of damage to the UAV as a result of exposure to atmospheric environmental 
factors. A computational experiment was carried out to determine the complexity of calculating the probability of UAV 
damage as a result of atmospheric forcing under various environmental conditions: “moderate conditions” at a wind 
speed of 3 m/s, precipitation intensity of 0.8 mm/h, and air temperature of 5 °C; “very heavy conditions” of the external 
environment at a wind speed of 12 m/s, precipitation intensity 3.5 mm/h, and air temperature of —6 °C. It was 
established that the use of triangular membership functions to calculate the probability of damage to the UAV as a result 
of atmospheric forcing provided higher performance compared to the rest ones (trapezoidal, pentagonal, Gaussian). 
Discussion and Conclusions. The values of the UAV damage probabilities obtained during the implementation of the 
methodology under atmospheric environmental influences can be used at the pre-flight preparation stage and during the 


flight to assess the feasibility of further performance of the flight task. Based on the analysis of the application of the 
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considered forms of membership functions, recommendations for their application are given. The use of algorithms with 


triangular membership functions will provide high performance of UAV control systems (UAV CS). 


Keywords: unmanned aerial vehicle, membership function, fuzzy sets, mathematical model, control, atmospheric 


forcing. 
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Introduction. The analysis of the experience in using UAV has shown that it is critical to increase the effectiveness 
of UAV under the conditions of environmental influences [1-4]. The equipment of modern unmanned aerial vehicles 
makes it possible to measure and assess the condition of technical systems, flight parameters, and the external 
environment [2]. 

With the rapid updating of current flight information and its large volume, the implementation of management faces 
the problem of providing the stable functioning of the UAV under the influence of the external environment. It is worth 
noting that human abilities are limited, and in a rapidly changing environment, a number of parameters of the UAV and 
the external environment, such as temperature and engine speed, wind speed and precipitation intensity, due to limited 
human capabilities, may not be taken into account or overlooked. Such control inaccuracies can cause the failure of the 
UAV equipment or its damage [3]. The influence of environmental factors on the UAV can significantly reduce the 
quality of the task, and for certain values of the parameters of external influences, completely exclude the possibility of 
using the UAV [4]. 

In this regard, solving problems related to maintaining the effective functioning of UAV under various atmospheric 
conditions is an important step in drawing up requirements for on-board UAV equipment that ensures flight safety in 
emergency situations arising under the influence of various weather phenomena [4]. 

Materials and Methods. The authors see a solution to the problem in automating the process of assessing the 
parameters of the external environment to determine the probability of damage to the UAV and in developing control 
corresponding to the current situation. In papers [1—5], meteorological conditions, specifically, precipitation, wind 
loads, temperature, air pressure and humidity, are considered as destabilizing effects of the external environment. To 
describe the atmospheric effects of the external environment', a probabilistic model is used in combination with the 
theory of random processes. A method is also known for estimating the intensity of precipitation and the water content 
of clouds based on remote measurement data’. The authors also note that the natural variability and stochastic nature of 
atmospheric conditions do not enable to establish strict mathematical dependences for describing atmospheric 
conditions of the environment. 

With the uncertainty of parameters and incomplete information about the state of the external environment, it is not 
always possible to present its unintended impact in the form of clear values. This problem can be solved using fuzzy set 
theory. In the case of its application, there is no need to know the exact mathematical model of atmospheric effects of 
the external environment to assess its impact on the UAV’. Assessment of the operating conditions of complex systems 


is one of the areas where fuzzy sets are widely used to eliminate the inaccurate nature of environmental influences [6]. 


' Kuznetcov IE, KuznetcovAA, Baklanov IO, et al. The Stochastic Model of Search and Detection of Ground Objects Using Unmanned Aerial 
Vehicles under Conditions of Irregular Influence of the Environment. In: Proc. VI Int. Conf. and Youth School “Information Technologies and 
Nanotechnologies (ITNT-2020)”. In 4 vol. VA Sobolev, ed. Samara: Samara National Research University; 2020. P. 152-156. 

* Kuznetsov DV, Bulgin DV. Metodika otsenki intensivnosti osadkov i vodnosti oblakov na osnove dannykh distantsionnykh izmerenii. In: Proc. XXI 
Int. Research-Methods.Conf. “Informatics: Problems, Methods, Technologies”, 11-12 Feb., 2021. Voronezh: Obshchestvo s ogranichennoi 
otvetstvennost'yu “VehlIborn”; 2021. P. 420-426. (In Russ.). 

> Cherepanov IE, Kuznetsov DV. Model’ prognozirovaniya stepeni slozhnosti meteouslovii s uchetom teorii nechetkoi logiki. In: Proc. XXI Int. 
Research-Methods.Conf. “Informatics: Problems, Methods, Technologies”, 11-12 Feb., 2021. Voronezh: Obshchestvo s_ ogranichennoi 
otvetstvennost'yu “Vehlborn”; 2021. P. 659-663. (In Russ.) 
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Thus, the analysis of well-known research in this area’. shows that in automated forecasting and decision support 
systems under conditions of lack and uncertainty of the initial information for its processing, fuzzy logic methods have 
proven themselves well. 

Research Results 

Development of methodology for determining the probability of damage to UAV as a result of exposure to 
atmospheric environmental factors 

In [6], a feature is noted that distinguishes the intelligent automatic control system from the one built according to 
the “customary” scheme. It is associated with the use of knowledge processing mechanisms to perform the required 
functions under uncertain (or incompletely specified) conditions with the random nature of external disturbances. 
Intellectualization of the UAV operation process to account for uncertainty factors can be achieved through the use of a 
knowledge base and an inference engine [6]. 

To determine the value of the probability of damage to UAV as a result of environmental factors, it is required: 

1. To formulate fuzzy estimates of atmospheric effects: 

— wind speed; 

— precipitation intensity; 

— temperature conditions of the environment; 

2. To develop a fuzzy model for calculating the probability of UAV damage as a result of atmospheric 
environmental impacts. 

3. To create a base of rules for fuzzy inference of the influence of atmospheric effects on UAV; 

4. To implement fuzzy inference. 

5. To defuzzify the fuzzy value of the probability of damage to UAV as a result of atmospheric influences of the 
environment Pap. 


The functional diagram of the technique is shown in Figure 1. 


1. Input of initial data 
(values of wind speed, precipitation intensity and air temperature) 


2. Assessment of UAV application environment conditions: 


2.1. Formation of fuzzy estimates of wind impacts: E\,(t) 


2.2. Formation of fuzzy estimates of atmospheric precipitation: G, (t) 
2.3. Formation of fuzzy estimates of the temperature conditions of the environment: Q,(t) 


3. Determination of probability of damage to UAV as a result of atmospheric influences: 


3.1. Development of a model for calculating the probability of UAV damage: [,(t) 
3.2. Development of an algorithm for determining the probability of UAV damage 


4. Output of results to the operator 


Fig. 1. Functional diagram of the methodology for determining the probability of damage to UAV as a result of 


exposure to atmospheric environmental factors 


*Ibid. 
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It is required to determine the probability of damage to UAV based on the analysis of the input fuzzy situation. Let 
us consider in detail the stages of the proposed methodology. 

Assessment of UAV application environment conditions 

Atmospheric environmental factors (wind, precipitation, and air temperature) may change during the flight. The 


dynamic change of the environment parameters has a significant impact on the stability of the UAV operation [7]. A 
short-range UAV made according to the aerodynamic scheme of an airplane type with flight speeds in the range from 
20 to 60 m/s is considered as a typical one. Icing of the surfaces of such UAV can cause a change in its mass, as well as 
aerodynamic characteristics. The wind speed can vary depending on altitude and flight time. It has a significant impact 
on the UAV of the class under consideration [8]. For this reason, it is reasonable to assess atmospheric impacts on the 
basis of dynamic fuzzy approaches. The initial data for the functional model are the environmental conditions: wind 
speed e, precipitation intensity g, and air temperature t. 


To estimate the effect of the wind acting on the UAV at time f, we introduce a set of triangular fuzzy numbers: 


E(t) = {Ex (), £20), 1 By), 1 Evga}, (1) 
where E(t) — a set of triangular fuzzy numbers describing the correspondence of the current wind speed e to fuzzy 
estimates of wind effects on UAV at time ¢; E(t) — a triangular fuzzy number describing the correspondence of the 
current wind speed e to a fuzzy estimate of wind impacts on UAV at time f; w = 1, New. 


In accordance with the concept of a fuzzy set, triangular fuzzy number E\,(t) can be understood as [9, 10]: 


Ey) = {(e, tg, (e))}, (2) 
where He (t) (e) — membership function of the wind speed e to the assessment of the wind influence E,,(t) on the 
UAV. 

Within the framework of the developed model, the qualitative assessment of the wind effect on the UAV depends 
only on its speed, i.e., possible inhomogeneities of the wind flow and its direction are not taken into account. 

The membership function of triangular fuzzy number uw B,(t) (e) can be represented by a tuple of three elements: 

a (©) = (le, (0): ¢e,(t),78,(O)), 3) 
where cg,,(¢) — clear value of wind speed that most accurately corresponds to the qualitative assessment E,,(t) of the 
wind influence on the UAV at time f; lg, (t) and rg,,(t) — clear values of wind speed in the least nonzero degree 
corresponding to the qualitative assessment &,,(t) of the wind influence on the UAV at time f. Values 
lz (t) and rg, (t) determine the fuzziness of estimate E.,(t). If le (t) = rg, (t) = cz, (t), then the estimate becomes 
clear. 

Specifically, £,(t) may mean an insignificant assessment of the wind impact on the UAV. If experts decide that the 
maximum wind speed at which its impact on the UAV under the current flight conditions at time ¢ can be interpreted as 
insignificant 0.8 m/s, then the fuzzy number membership function will take the form He ct) (e) = (0; 0; 0.8). Since the 
first and second elements of the tuple are equal to zero, the region of nonzero membership is located on a closed 


interval [0, 0.8]. To determine the degree of membership He ct) (e) of the wind effect at a speed of e to the fuzzy 


assessment of the impact E£,(t), it is required to substitute e into the equation of a straight line passing through the 
points (0;1) and (0.8;0). 

During the flight, UAV can determine the actual wind speed offline. Based on the values of tuples of three elements 
set by experts, by analogy with the above example, we formulate possible estimates of the influence of wind on the 
UAV. For this, the values of the triangular membership functions of the fuzzy sets F(t) from the values of the wind 


speed e are represented by graphs (Fig. 2): 
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— E,(t) — “insignificant influence”, He (2) = (0; 0; 0.8); 
— E(t) — “very weak influence”, He cr) (€) = (0; 1;3); 

— E(t) — “weak influence”, HE (t) (e) = (0.2; 2.5; 4.8); 

— E(t) — “strong influence”, He ct) (e) = (1.5; 4.2; 6.9); 

— E(t) — “very strong influence”, Het) (e) = (1.5; 7; 12.5); 
— E(t) — “dangerous influence”, HE ct) (e) = (5; 11; 18); 


— E,(t) — “very dangerous influence”, HEC) (e) = (10; 15; 00), 


V\/ 
y 
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Fig. 2. Graphs of dependences of triangular membership functions of fuzzy 
sets Ey, (t) on the values of wind speed e 


By analogy with the assessment of the influence of wind, we consider estimates of the influence of precipitation and 
temperature affecting the UAV at time f. 


To assess the influence of precipitation, we introduce a set of triangular fuzzy numbers: 
G(t) = {G,(0), G00), Gn (O), ns Green (t}, (4) 
where G(t) — a set of triangular fuzzy numbers describing the correspondence of the current precipitation intensity g 
to fuzzy estimates of the effects of precipitation on UAVs at time f, G,(t) — triangular fuzzy number describing the 


correspondence of the current precipitation intensity g to a fuzzy estimate of the effects of precipitation on UAV at time 


t h= 1, Nee): 


To assess the effect of temperature, we introduce a set of triangular fuzzy numbers: 
Q(t) = {0,0 G20), 1 O(0), Ono (OF, (5) 
where Q(t) — a set of triangular fuzzy numbers describing the correspondence of the current temperature q to fuzzy 


estimates of the effects of temperature on the UAV at time ¢; Q,(t) — triangular fuzzy number describing the 


Mechanics 


197 


http://vestnik-donstu.ru 


198 


Advanced Engineering Research 2022. V. 22, no. 3. P. 193-203. ISSN 2687-1653 


correspondence of the current temperature qg to a fuzzy estimate of the temperature effects on the UAV at time 4; 
s= Tow: 

In the developed model, the influence of weather conditions is estimated based on three sets (t), G(t), and Q(t). 
The elements of these sets correspond to triangular membership functions, which simplifies the adaptation of the model. 
To adjust the membership function, it is only required to change the values of atmospheric exposure that correspond to 
the greatest and/or least extent to the adaptable fuzzy estimate. 

Algorithm for determining the probability of UAV damage as a result of exposure to atmospheric 
environmental factors 

To implement the methodology, an algorithm has been developed to determine the probability of damage to UAV as 
a result of exposure to atmospheric environmental factors. The algorithm enables, under conditions of uncertainty, to 
take into account the influence of atmospheric environmental factors to calculate the probability of damage to UAV 
according to the current situation. The initial data for the algorithm are fuzzy sets E\,(t), G,(t) and Q,(t). 


Figure 3 shows its block diagram. 
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Aggregation of the Output of results 
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Fig. 3. Block diagram of the algorithm for determining the probability of damage 
to UAV as a result of exposure to atmospheric environmental factors 


Determination of the probability of UAV damage as a result of atmospheric environmental influences 
Papers [4, 5] note the natural variability and stochastic nature of atmospheric conditions. Due to the fact that the 


atmospheric effects of the external environment are dynamic, it is reasonable to take this circumstance into account 
when constructing a model for determining the probability of damage (equipment failure) of UAV. We introduce a set 


of triangular fuzzy numbers: 


10) = {0,40 £0, tj Of, (6) 
where I(t) — a set of triangular fuzzy numbers describing the correspondence of the current value of the probability of 


damage to UAV Pap at time r; [,(¢) — triangular fuzzy number describing the correspondence of the current value of 
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the probability of damage to the UAV Pag to a fuzzy assessment of the level of exposure to atmospheric environmental 
factors on UAV at time f; z = 1, Nia): 

With regard to the UAV performance levels [11] (operable, limited operable, and inoperable condition), consider the 
following UAV damage probabilities: [,(t) — “low” (approximately 0.2), [,(t) — “medium” (approximately 0.5) and 
[,(t) — “high” (approximately 0.8). Triangular fuzzy number [,(t) is represented as [9, 10]: 


,(t) = {(Pap» Hz, (Pas) }.2 = 13, (7) 
where Le (ty (Pap) — the membership function of the probability of damage to UAV Pg to the assessment of impact 


7, (t) of the wind on the UAV. 

Given that the model under consideration has a dynamic character, the compilation of the fuzzy inference rule base 
should be preceded by a study of the influence of atmospheric conditions on the probability of damage to UAV. To that 
end, a series of experiments were conducted under various environmental conditions. Based on expert assessments of 
statistical data, a database of fuzzy inference rules was formed [10]. The rules for determining qualitative estimates of 
the probability of damage to UAV as a result of atmospheric influences of the external environment were formalized 
using fuzzy sets as follows: 

If gis Q,(t) and eis £,,(t) and gis G,,(t) then is /,(t). (8) 

That is, under simultaneous conditions of minimally sufficient membership: 

— of temperature g to the fuzzy assessment of the impact of temperature Q,(t); 

— of wind speed e to the fuzzy assessment of the wind impact E,, (t); 

— of precipitation intensity g to the fuzzy assessment of the impact of precipitation G, (t), the fuzzy assessment of 
probability L(t) becomes relevant. 

To obtain clear value Px,(t) of the probability of damage to UAV as a result of atmospheric influences of the 
external environment at time ¢, it is required to defuzzify the corresponding fuzzy result according to formula [10]: 

f Pap rar (Pap)d Pap 


Pxp(t) = 1 a = 
0 Mee) (Pap)@Pap 


(9) 


3 


where Me (Pap) — the resulting membership function of fuzzy set [,, representing the general conclusion (conclusion) 


from all the rules (8); ,, — the output variable of the fuzzy inference, characterizing the fuzzy value of the probability 
of damage to the UAV as a result of atmospheric environmental influences. 

Thus, the intellectualization of the functioning process is achieved through the use of a knowledge base and an 
inference engine, and provides that uncertainties are taken into account [6]. 

Scheme of the UAV operation process under environmental conditions 

To calculate the probability of damage to UAV as a result of atmospheric environmental influences, it is required to 
implement fuzzy inference. It consists in the sequential execution of aggregation, activation and accumulation 


operations [10]. The scheme of the UAV operation process under environmental conditions is shown in Figure 4. 
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Algorithm for determining the probability of UAV damage as a result of exposure to atmospheric environmental factors 
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Fig. 4. Scheme of the UAV operation process under the influence of atmospheric environmental factors 

Computational experiment 

In [12], an analysis of foreign periodical literature in recent years has been carried out to identify formal approaches 
to dealing with uncertainties. The conclusion was made about the expediency of using n-angle numbers to determine the 
accuracy of the uncertainty formalization. Thus, when using the triangular representation of fuzzy numbers [10], each 
initial fuzzy number is described by three scalar values, which greatly simplifies the computational process. 

In many cases, triangular and trapezoidal representations of fuzzy numbers are insufficient. If it is required to 
consider more complex influence of uncertainties, specifically, due to the greater exposure of the equipment installed on 
UAV to precipitation, it is reasonable to use more complex forms of membership functions [12, 13]. Therefore, in [13], 
pentagonal fuzzy numbers, which reflect uncertainties much more accurately, were studied as an alternative. 

Let us consider the basic, most used membership functions for the representation of fuzzy linguistic terms: 
Gaussian, triangular, trapezoidal [10], and pentagonal fuzzy numbers [13, 14]. 

Based on the analysis of papers [13, 14], we present the pentagonal membership function of elements to fuzzy sets 


of wind load effects E,,(t) in the form of a tuple: 

He Ce) = (le, 3 ke, (0; ca, (O; pe, (O72, 0), (10) 
where cg, (€) — clear value of the wind speed most accurately corresponding to the qualitative assessment E,,(t) of the 
wind impact on UAV at time #; lg (t) and rg, (¢) — clear values of wind speed in the smallest (nonzero) degree 
corresponding to qualitative assessment E,,,(t) of the wind impact on UAV at time 1; kz, (t) and pg, (t) — clear wind 
speed values with a degree (nonzero), corresponding to qualitative assessment F(t) of the wind impact on UAV at 
time ¢. Values kg, (t) and pg, (t) determine the fuzziness of assessment E(t). In case, if le (t) = rg, (t) = 
kz, (t) = pz,,(t), we get a triangular fuzzy number. By changing the values of « level, we obtain a different form of 


the pentagonal membership function (Fig. 5). 
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Fig. 5. Pentagonal membership function in MATLAB software environment 


To study the pentagonal membership function in a software environment, it is required to describe it in the form of a 
program code. When implementing the algorithms*”, a custom pentagonal membership function was developed in the 
Fuzzy Logic Toolbox application package of the MATLAB software environment. This tool provides exploring 
mathematical models and algorithms based on the application of fuzzy logic theory. 

To assess the accuracy of uncertainty when working with the described membership functions, it is required to 
conduct a computational experiment. It is reasonable to carry out the comparison on the same set of source data. The 
values of the probability of UAV damage as a result of atmospheric influences are calculated under “moderate 
conditions” (actual values: wind speed 3 m/s, precipitation intensity 0.8 mm/h, air temperature 5 °C) and “very difficult 
conditions” (actual values: wind speed 12 m/s, precipitation intensity 3.5 mm/h, air temperature — 6 °C) of the external 
environment. 

Using MATLAB tools, the performance of the software package implementing the developed algorithm* was 
evaluated depending on the form of the membership function used. An J/BM-compatible PC with minimal system 
requirements (P4 — 2.500 MHz, 2.048 Mb RAM, Windows operating system) was used to conduct the computational 
experiment. Triangular, trapezoidal, pentagonal and Gaussian membership functions were evaluated. The results are 
presented in Table 1. 


Table | 
Average operating time of the software package 
Time, s 
Type of membership “Moderate conditions”: “Very difficult conditions”: wind speed 
function used wind speed 3 m/s, precipitation intensity 12 m/s, precipitation intensity 3.5 mm/h, 
0.8 mm/h, air temperature 5 °C air temperature — 6 °C 
Triangular 0.028 0.028 
Trapezoidal 0.031 0.03 
Pentagonal 0.034 0.035 
Gaussian 0.031 0.03 


From Table 1, it can be concluded that a computer program based on the use of triangular membership functions has 
a higher speed in comparison to the others. 

A promising area of research is studying the influence of the form of the membership function on the accuracy of 
calculations of the probability of damage to UAV under the uncertainty of the impact of atmospheric factors of the 
environment on the functioning of UAV. 


* Belonozhko DG, Korolev ID, Chernyshev Yu.O. Certificate of State Registration of Computer Program No. 2022613419, Russian Federation. 
Software Package for Evaluating the Accuracy of Uncertainty and Performance of Fuzzy Algorithms for Controlling an Unmanned Aerial Vehicle, 
no. 2022613002; 2022. (In Russ.) 

° Belonozhko DG, Korolev ID, Chernyshev Yu.O. Certificate of State Registration of Computer Program No. 2022613691, Russian Federation. 
Software Package for Intelligent Control of an Unmanned Aerial Vehicle under Conditions of Destructive Impacts Based on Pentagonal Accessory 
Functions, no. 2022613009; 2022. (In Russ.) 
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Discussion and Conclusions. A method has been developed for determining the probability of damage to UAV as a 
result of exposure to atmospheric environmental factors, based on the mathematical apparatus of fuzzy sets. Fuzzy 
qualitative estimates made it possible to significantly expand traditional mathematical modeling techniques that require 
accurate information about input quantities. The use of this technique provided assessing the probability of UAV 
damage when the outcomes of alternatives were not known accurately and the probabilities of their occurrence were 
estimated using membership functions. The obtained values of the probabilities of UAV damage as a result of 
atmospheric influences of the external environment can be used at the stage of pre-flight preparation and during the 
flight to assess the feasibility of further performance of the flight task. 

With the help of a software package implementing this technique, the values of the probability of damage to UAV 
under “moderate conditions” and “very difficult conditions” of the external environment were obtained. According to 
the results of the analysis of the application of these forms of membership functions and the computational experiment, 
the high performance of the software package using the triangular method of formalization of fuzzy sets was validated. 
A significant advantage of triangular membership functions is that their definition requires the least amount of 
information compared to other functions. 

Thus, it is reasonable to apply algorithms using triangular membership functions in the onboard computers of the 
UAV CS to assess the atmospheric parameters of the environment. This will provide high performance of the UAV CS 
under conditions of using computers with low performance. 
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Abstract 
Introduction. The problem of controlling the lower-extremity powered exoskeleton motion was investigated. To solve 


it, it was proposed to use a program control and feedback control. The formation of control in the form of feedback 
required an assessment of the state of the exoskeleton (rotation angles, angular velocities, and accelerations of the 
links). The possibility of using an inertial measuring unit to estimate angular velocities and accelerations of exoskeleton 
links was considered. The work objective was to develop laws for the formation of the exoskeleton motion control, 
which could provide the stability of the program motion and use the measurements of encoders, micromechanical 
gyroscopes and accelerometers. 

Materials and Methods. Previously performed mathematical modeling of the exoskeleton dynamics was used to form a 
program control. It was proposed to equip the exoskeleton with inertial sensor units. This solution made it possible to 
evaluate the state vector of the exoskeleton and to use these estimates in a feedback loop. A mathematical model of 
measurements of these sensors was described. The proposed version is suitable for control systems of three-link 
exoskeletons of the lower extremities and can be expanded to the case of multi-link exoskeleton designs. 

Results. New laws of exoskeleton motion control based on a mathematical model of the system dynamics and using 
measurement information from encoders and inertial information sensors were proposed. Numerical simulation of 
exoskeleton motion was performed in the Wolfram Mathematica mathematical package. Its results confirmed the 
operability of the proposed control and the possibility of using an inertial sensor unit to assess the exoskeleton state. 
The numerical simulation results for the following program movements were presented: lifting the exoskeleton from a 
sitting position to a vertical position, and stabilization of the vertical equilibrium position. 

Discussion and Conclusions. The proposed control can be applied in exoskeletons for medical purposes, e.g., in the 
task of verticalization of patients with dysfunctions of the musculoskeletal system. The possibility of using 
measurement information obtained from inertial measurements units in the problem of estimating the state of 
exoskeleton links was demonstrated. The use of inertial sensors will make it possible to determine the angular 
acceleration of the exoskeleton links, avoiding numerical differentiation of the measurement information received from 
the encoders. The obtained estimates of angular acceleration allow us to introduce feedback on angular accelerations 


into the control system, which opens up the possibility of improving transients in controlling the exoskeleton motion. 
Keywords: exoskeleton, mathematical model, inertial sensors, control, state estimation. 
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Introduction. The development of devices designed to facilitate motion and increase motor efficiency (of 
exoskeletons) is perspectives for development robotics. 

Exoskeletons are increasingly used in industry [1], military [2], medicine [3], and other fields. This specifies 
topicality of the research on their circuit and structural design, energy consumption optimization, motion trajectory, 
dynamics, control. 

The present paper considers an active lower-extremity powered exoskeleton of a human with a rigid structure of the 
load-bearing frame. It is used for the rehabilitation of the musculoskeletal patients. 

The development of such exoskeletons requires solving interrelated scientific and technical problems. They include: 

— formation of a skeletal scheme; 

— mathematical modeling of exoskeleton motion [4-6]; 

— optimization of the links motion based on a mathematical model of kinematics and dynamics’; 

— reasonable construction of the motor system [7, 8]; 

— reasonable construction of the motor system [9—11] and estimation [13, 14]. 

All of the above tasks are connected by the problem of finding control actions for program motion. Firstly, their 
source may be the dynamic model of the system. Secondly, they are developed in the form of a proportional integral- 
differential (PID) controller used to stabilize the program motion. A combination of two methods is also possible. At 
that, for the formation of feedback control (in particular, in the form of a PID controller), algorithms for assessing the 
state of a dynamic system and a set of sensors whose indicators are used in algorithms, play an important role. 

To obtain information about the position of the exoskeleton links, they are equipped with an encoder or gyroscope 
[9-14]. Other possible options include accelerometer, magnetometer, potentiometer, sensors of torque, force, biosignals, 
etc. 

In [13], it is proposed to use one block of inertial sensors consisting of gyroscopes and accelerometers on each link 
of the exoskeleton. This solution will provide estimating the angles of rotation, angular rates and accelerations. At the 
same time, the measurement model turns out to be nonlinear, which complicates the structure of control and monitoring 
systems. 

The study aimed at improving the control of the exoskeleton by using two inertial blocks attached to each link of the 
exoskeleton. This approach enables to estimate the rotation angles, angular rates, and accelerations of the exoskeleton. 
With the help of the proposed arrangement of gyroscopes and accelerometers, a linear measurement model can be 


obtained, which will open up ways to simplify the control system and apply the theory of optimal evaluation of linear 


' Yatsun SF, Savin SI, Yatsun AS, et al. Optimization of Time of Verticalization of the Exoskeleton by Criterion of Energy Efficiency. In: 
Vibratsionnye tekhnologii, mekhatronika i upravlyaemye mashiny. Pt. 2. Kursk: Izd-vo Yugo-Zap. gos. un-ta; 2016. P. 151-160. (In Russ.) 
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systems [15]. This, in turn, will be the basis for improving the motion control of the exoskeleton. The criterion for 
improving the control quality is reducing the deviations of the angular rates of the links from the program movement. 
Materials and Methods. Consider the kinematic scheme of an exoskeleton with a rigid structure of the power frame 


with symmetrical movement of the legs in the sagittal plane (Fig. 1). 


Fig. 1. Kinematic diagram of the lower-extremity powered exoskeleton in the sagittal plane 
C2, C3, Cs — centers of mass of the lower legs, thighs and body, respectively; A), 42, 43 — joints connecting the 
links of the exoskeleton; @2, 3, @4 — rotation angles of the exoskeleton links; M2, M3, M, — control moments created 
by the actuators located in the joints A), A2, A3, respectively. Coordinate system xyz is introduced to describe the motion. 
x-axis is directed along the reference surface, y-axis is directed along the local vertical to the reference surface. z-axis 
complements the coordinate axes to the right Cartesian coordinate system. 


As a basis for building an exoskeleton motion control system, we use a dynamic model from [6]. 


A(q)4+F(q)q° +Dq+P(q)=M. (1) 
Here: 
Jy J, cos(@, —9;) J, C0S(, -@,) 

A(q)=| J,; cos(@, —9;) le J,, cos(@;—@, ) | — inertial force matrix; 

J,, cos(@; —@,) Jz, C08 (9; —9,) Jag 

0 J,, sin(@, —@;) J,,sin(@, —,) 

F(q) =| —J,, sin(~, —9;) 0 J,, sin(@; —,) | — speed force matrix; 

Soy sin(9, ~9,) Soy sin(@, —9,) 0 
D= diag(u,, L,5[ a) — diagonal matrix of dissipative forces characterizing linear friction in joints; 
P(q) = (M,. cos~, M,,cos@, M,, cos, y' — column vector of moments arising from the action of gravity; 


M= (M, —-M, M,-M, M, y — column vector of control torques generated by actuators located in the joints of 


Saypulaey G. R., et al. Application of an Inertial Sensor Unit for Position Estimation and Motion Control of the Lower-Extremity Powered 


the exoskeleton links; q = (9, QP, Dy i — column vector of generalized coordinates, consisting of rotation angles of 
exoskeleton links; J jk ( ik =2, 3,4) — moments of inertia of exoskeleton links; 
Hy (k — 2,3,4) — coefficients of linear friction in the joints of the exoskeleton links; M,, (k = 2, 3,4) — maximum 


values of the moments of gravitational forces perceived by the output links of the drives. 


Let the laws of change of angles and their derivatives be given for program motion: 
PF =O, GF =O, G7 =). (1=2,3,4). (2) 
To control the exoskeleton motion, it is proposed to create control actions in the form of a sum of control moments: 


— with software management M/’ = M, (9° 67.67 ) (i = 2,3, 4) : 


—PID controller M;"” =M,(@/ —9,,6/ —6,.6 4, ),(i=2.3,4). 
M, =M3+M3", M,=M?+Mi”, M,=M;+M)”. (3) 
For program control, the values of control actions M/ can be calculated using the formulas [7]: 
Mi =Mg. cose +194 +5.G! +Jss| 8 cos(o% -o%)-(98)' sin(o! -0%) |+ 
+Jua| Gf c0s( of -o1)-(62) sin(o’ 0) | 


M3 =Mg; Coss +1565 + 5,65 + Js E cos(9} — 9; )-(@3) sin(e} -9)|+ 


(4) 
“i .p\2_. 
+ Jue] GF c0s( of -0!)+(6) sin(of -9') |+ Mi, 
Me=M P aes eee a Pop »p\? - ( pop 
2 = Mg. COS, FHP, +I,Q. +453) Pz COS(P, —P3 J +(P; } SIN|P, —; } |+ 
. <p\? 
+Jui| 6! c0s(! -of) +(0%) sin(o! ~9) | +M%. 
Values M/" are calculated according to the feedback principle (in the form of a PID controller): 
M’P =K f(g" b,)+Kp,(@) —@,)+K [(o2t,)-¢ (t,))dt 
4 D4 at D4 — D4 pa\P4 —P4 14 J (Pa) — Pah b 
0 
dyip_, Pp r/. : 
M;” =K,, dt (6; $s) Kp, (6; 6s) Ky, { (68 (t,) G3 (t,)) dt, + M2, (5) 
0 
dyip sp og Te 2 
M” =K,, races $,) Ky (65 -,)+K,, [(@5(t,) -@,(t,)) dt, + M2”, 
0 


or 


M;” =K,, (6: —6,)+ Koa (6; ~@,) Ky, (0 0,4), 
M;" =Kp, (6: —é,)+ Kp; (o; $3) +K 3 (9: 03) tM,” (6) 
M; =K,, (6; —G, ) +K,, (6; $2) Kp (9; 2) cM," 


Here, the coefficients of the PID controller Kp;, K;,, Kp; (i = 2, 3, 4) can be: 
— constant (they are found from the stability conditions); 
— time functions (found from solving optimal control problems [9, 13]. 
One of the problems of creation of control actions in the form of a PID controller is possible errors in the 


differentiation of measurement information. To avoid them, we use: 
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— encoders (for measuring the current rotation angles @, of the exoskeleton links); 
— inertial sensor blocks of micromechanical gyroscopes (for measuring angular rates @, ) and accelerometers (for 
estimating angular rates @, and accelerations 4, ). 
The model for measuring rotation angles using encoders can be presented as: 
PE2 =P2, Pes t+ Per =P3, Peg + Pe3 + Pps = Py, (7) 
where ,5,0 35,4 — indications of encoders located in joints A;, Az, A3 (Fig. 1). 
The model of angular rate measurements using micromechanical gyroscopes can be written as: 
Q67 =O, Q63 = 93, Qe =O, (8) 
where Q¢9, N63, Qe, — indications of gyroscopes attached to the links of the exoskeleton. 
To estimate angular rate and accelerations using accelerometers, a pair of two-axis accelerometers can be used, 


which are located at opposite ends of each of the links. They measure the apparent accelerations associated with the 
corresponding links in projections on x; y,-axis (i = 1, 2, 3). 
A A A A A A ‘ 2 a: 
Pala, =A = (ras, = as (i. avee ), = (103 13) : 
A A A A A A : 2 x 
we TR A = (Ges avis : Gen aye i. ~ (105 -165) 3 (9) 
c A c A c A : 2 r 
fac —fxvc = (Ge —fxvc ). ge —fxvc i = (1,0; -1,6,) : 
Here, fe, , fe. — vectors of readings of two-axis accelerometers mounted on link A,A, and located at points A, and 
"2 2 
A,, respectively; f; 4a) ie — vectors of readings of two-axis accelerometers mounted on link A,A, and located at 


points A, and A,, respectively; f aos ae — vectors of readings of two-axis accelerometers mounted on link A,C and 


located at points A, and C, respectively; J, =|A,A,|,/, =|A,A, 


,1,=|A,C| — distances between pairs of 


accelerometers attached to links A,A,, A,A, and A,C, respectively. 


We combine equations (7)-(9) and write formulas for estimating the angles of rotation, angular rates, and 


accelerations: 
> =Per> 93 =Pest+ Ger, Ps =PpatPe3+ Op, G2 =QG, G3 =QG3, Oy =QGu, 
(te,-the), (Bef), (Feet ao 
1 os 2 6 


a ¥3 
9, 1 » Qs L » D4 1, : 


To correct estimates of angular speed modules, the equations can be used: 


(Ey ~ Lara ). = 1,63, (i 7 re ), = L465, (fe a ie) - 1,04. (11) 


x3 
Thus, to estimate the rotation angles, angular rates and angular accelerations of the exoskeleton links, measurements 
can be processed with a moving average filter, Kalman filter algorithm, etc. [10, 13, 15]. 
In the simplest version of closing the feedback circuit when forming the control, we obtain expressions for the 


control moments: 
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(fc fat) 


M,=M;,+K,,| 4; 4 = 1+Kpy (6 05) Ki, (9: rs — Pes ®x2)> 


I, 
(ot 
M, =M; +K,,| 6; 4 i = |+Kos (6; Os) KK, (95 9x3 2) +My”, (12) 
2 
fh fA 
Mea + gb +O toh + Kya ($ Oca) Kia (0! 2) + ME. 
1 


The introduction of angular acceleration feedback into the control system provides the diagonalization of the 
inertial forces matrix to reduce the cross-influence between generalized coordinates. 

We simulate the exoskeleton motions to test the operability of the proposed control law. As a program motion, 
consider: 

— change of position from sitting to standing [7]; 

— stabilization of the vertical position of the exoskeleton with initial deviations of the links from the vertical. 

Research Results. For numerical simulation of the system, we took the parameters of the mathematical model of 
the exoskeleton from [6, 7]. We considered a special case by choosing the same coefficients of the PID controller for 
each drive: K,, =1,000N-m-s, K, =1,000N-m, K,,=0N-m/s(i=2,3,4). Figure 2 shows the time dependences 
for the rotation angles and angular rates of the exoskeleton links. They were obtained at control moments corresponding 
to rising from a sitting position. As can be seen from Figure 2 a, at the end of the motion, the rotation angles of each 


link are equal ©, =@, =@, =90° , which corresponds to standing. 


Q;° 


150 


100 


50- 


a) b) 
Fig. 2. Simulation results of rising from a sitting position without initial deviations from the program motion in the angles 


of rotation of the exoskeleton links: a — rotation angles of the links; b — angular rates of the links 


Now let us consider the case with deviations of about 20° in the rotation angles of the links compared to the initial 
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position set in the program motion. The simulation results of rising from a sitting position are shown in Figure 3. As in 
Figure 2, at the end of the motion time, the values of the rotation angles were equal to 9, =@,=@, =90°. 


Consequently, the use of a combination of software control and control in the form of a PID controller accurately 


brought the exoskeleton to the desired position. 209 
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6") 


150) 
100} 


50} 


a) b) 
Fig. 3. Simulation results of rising from a sitting position with initial deviations from the program motion in the rotation 
angles of the exoskeleton links: a — rotation angles of the links; b — angular rates of the links 
Thus, when modeling program motion, the operability of the proposed law of formation of control actions described 


in (4) and (12) was shown. 


The vertical position of the exoskeleton is unstable without control; therefore, disturbing factors can cause 


deviations from equilibrium up to the falling of the exoskeleton. 


We checked the operability of the proposed control in the task of stabilizing the vertical position of the exoskeleton. 
Values for the program motion were: 9% (t) = @4 (t) = 94 (t) = 90° . We considered the case when there were deviations 
of about 20° from the vertical position in the rotation angles of the links (Fig. 4). 


9;.° 
120+ 


o/s 


100} 


a) b) 
Fig. 4. Simulation results of stabilization of the vertical position with initial deviations from the program motion in the 


rotation angles of the exoskeleton links: a — angles of rotation of the links; b — angular rates of the links 
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Figure 4 shows that the proposed control provides stabilization of the vertical position of the exoskeleton links. 


For comparison, we present the results of modeling the stabilization of the vertical position using the following 


values of the coefficients of the PID controller: K,,=1,000 N-m-s, K, =1,000N-m, K,, =100 N-m/s (i= 2,3,4) 
(Fig. 5). 

@,.° 

120 


100: 


80} 


60! 
40; 


20 


a) b) 
Fig. 5. Simulation results of stabilization of the vertical position with initial deviations from the program motion in the rotation 


angles of the exoskeleton links: a— angles of rotation of the links; b — angular rates of the links 


It can be seen from Figure 5 that when using the differential link ( K,, #0) of the regulator, it was possible to 


reduce emissions in the dependences of the angular rates of the exoskeleton links compared to the case K,, =0 
(Fig. 4 5). 

Thus, the possibility of controlling the lower-extremity powered exoskeleton links using the readings of encoders 
and blocks of inertial sensors (micromechanical gyroscopes and accelerometers) was considered. The simulation results 
have validated this approach. 

Discussion and Conclusions. The proposed solution to the problem of controlling the motion of the exoskeleton 
can be used in medical applications, including for the verticalization of musculoskeletal patients. 

Thus, the data of the inertial sensor blocks can be used in the task of assessing the state of the exoskeleton links — 
their angular rates and accelerations. At the same time, there is no need for numerical differentiation of the measuring 
information of the encoders. Therefore, feedback on angular accelerations can be introduced into the control system, 


which makes it possible to improve transients when controlling the exoskeleton motion. 
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Abstract 
Introduction. The problem of critical loads of a compressed orthotropic rectangular plate on an elastic base was 


considered. The following orthotropy parameters were set for the plate: Poisson coefficients, Young's modules for the 
main directions, and the shear modulus of the plate material. The components of the compressive load were uniformly 
distributed along two opposite edges of the plate and acted parallel to the coordinate axes. The edges of the plate were 
loosely pinched or pivotally supported. Cases were also considered when two parallel edges of the plate were free from 
loads, and the other two were freely pinched or pivotally supported. 

Materials and Methods. The problem was studied on the basis of a system of nonlinear Karman-type equilibrium 
equations. The critical values of the load parameter were determined from a linearized problem based on a trivial 
solution. At the same time, the variational method in combination with the finite difference method was used to solve 
the boundary eigenvalue problem. 

Results. The problem was reduced to solving a parametric linear boundary eigenvalue problem. In case of boundary 
conditions of a movable hinge support, exact formulas of eigenvalues and eigenfunctions were given. While in case of 
free edge pinching, a variational method was used in combination with a finite-difference method, and a computer 
program for solving the problem was built. It was established that one or two eigenfunctions expressing the deflection 
of the plate could correspond to the critical value of the compressive load parameter at which the stability of the 
compressed plate was lost. The results of numerical calculations of the critical values of the compressive load at 
different values of the orthotropy parameters were presented, and graphs of the corresponding equilibrium forms were 
constructed. For the case of a long orthotropic plate on an elastic base, it was established that the main term of the 
asymptotic expansion of the solution to the linear eigenvalue problem was determined from the problem of critical loads 
of a compressed beam on an elastic base with an elastic modulus that coincides with the elastic modulus of the plate in 
the longitudinal direction. 

Discussion and Conclusions. The problem of critical loads of an orthotropic plate compressed in two directions lying 
on an elastic base was investigated. As the compressive load component increased along one direction, the critical value 
of the load compressing the plate along the other direction decreased. If an orthotropic plate was compressed by a load 
along a direction that corresponded to a greater bending stiffness, then the critical value of the loss of stability was 
greater than the critical value of the compressive load acting along the direction of a lesser bending stiffness. The 


presence of an elastic foundation increased the bearing capacity of the compressed plate. 
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Introduction. L. M. Zubov in [1] derives the equilibrium equations of an elastic plate containing sources of internal 
stresses in the form of dislocations and disclinations, which are a modification of the Karman equations. In this work, 
the problem of bending a thin plate (membrane) under the action of internal stresses caused by defects is also solved. 
In [2-5], the theory of dislocations and disclinations is used in the study of the behavior of plates and shells. In 
article [6], a general solution of the oscillation equation of a rectangular orthotropic plate with free edges is constructed. 
In [7, 8], asymptotic solutions of problems of statics and dynamics of narrow plates are constructed. 

In [9], the problem of equilibria branching under the influence of a small normal pressure of a rectangular plate 
compressed in one direction with dislocations and disclinations is considered. The case is considered when the 
longitudinal edges of the plate are free from loads, and the other two edges are fixed. In [10], an asymptotic solution is 
provided for the problem of critical loads of a compressed narrow plate using a small parameter expansion of the 
relative width of the plate. In [11], the branching problem of the equilibria of an elastic rectangular plate with internal 
stress sources is investigated for the case when compressive loads are unevenly distributed along the edges. It is 
established that with even forms of the incompatibility function and even forms of distribution along the edges of 
compressive forces, the presence of a small normal load does not reduce the bearing capacity of the plate. 

In work of Morozov N. F., Belyaev A. K., Tovstik P. E., Tovstik T. P. [12], an asymptotic derivation of two- 
dimensional equilibrium equations of a thin elastic inhomogeneous plate made of a general anisotropic material is 
presented. Conclusions about the forms of loss of stability of a compressed plate on an elastic soft base are formulated 
in [13]. In [14], the bending stiffness of a thin elastic multilayer plate with transversely isotropic layers is investigated. 
Two models of accounting for the effect of transverse shear are compared. 

In [15], a model of thermoelastic viscoplastic deformation of a composite that is cross-reinforced with continuous 
fibers in arbitrary directions is constructed. 

In [16], within the framework of a direct approach to plate theory, the problem of a hyperelastic plate with 
inhomogeneously distributed initial stresses is considered. The plate is considered as a material surface with five 
degrees of freedom (three displacements and two rotations). In [17], the stability of a uniformly compressed circular 
two-layer plate with an initially compressed or stretched layer is analyzed within the framework of nonlinear elasticity. 
A model of incompressible neo-Hookean material is used for determining the material ratio. The analysis of the 
dependence of the resulting critical stresses on the initial deformations and stiffness parameters is presented. In [18], the 
bending of a three-layer plate with surface and interfacial stresses is investigated, the theory of plates with first-order 
shear deformations and the Gurtin-Murdoch model of surface stresses are considered. The dependence of the plate 
stiffness parameters on the surface elasticity modules is analyzed. In [19], the problem of instability of a three-layer 
nonlinear elastic rectangular plate with a prestressed middle layer is considered. 

In this paper, the influence of the parameters of the plate material and the elastic base parameter on the critical loads 


of the loss of stability of a compressed elastic rectangular orthotropic plate on an elastic base is investigated. 
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Materials and Methods. Let us consider an orthotropic rectangular plate lying on a linearly elastic base, which is 


compressed by forces P and Q along the X and Y axes. The equilibrium equations can be written as [7]: 


D,0,. W + 2D,0;,.0, W + D,0,W+KW =[W,F], 


2 (1) 
 otpy| £—-™ leeatrs + atr=—[w,w]. 
E, G E£, E, 2 
Let us write the boundary conditions in the form: 
W=0,W=0, 6,0,F=0, 60;F=-Pat |X|=a/2, bs 
W=0,W=0, 6,0,F=0, 0,F=-Qat |¥|=b/2. 
W=6,W=0, 6,0,F=0, 0,F=-Pat |X|=a/2, és 
W=0,W=0, 0,0,F=0, 0:F=-Qat |Y|=b/2. 
W=0,W=0, 0,0,F=0, O0;F=—P at |X|=a/2, 4) 
0, W +v,0;,W = 0, W +(2-v,)0,0,W =0,0,F = 6, F =0 at |Y|=b/2. 
W=6,W=0, 0,0,F=0, O0;F=—-P at |X|=a/2, 6) 
0, W +v,0;,W = 0, W +(2-v,)0,0,W =0,0,F = 6,F =0 at |Y|=b/2. 


[W,F] = 05, W0,F+ 0, W0;,F-26,0,W0,0,F, D,=Eh’/(12(1-v,v,)), i=12; D, =Gh’/12, D,=D,v,+2D,, 
E,v, =E,y, . 


Here v,, v,, E,, E,, G— Poisson's coefficients, Young's modules for the main directions and the shear modulus 
of the plate material; D,, D, — bending stiffness in the main directions; D, — torsion stiffness; the function 
W(X, Y) expresses deflection, and F(X, Y) denotes the stress function; 4 — plate thickness, KW — base reaction; 


K — modulus of subgrade reaction; the origin of coordinates X, Y is in the center of the plate and the axes are parallel to 
its edges; 0 — partial derivative of the n-th order with respect to the variable X. We assume that the edges of the plate 


are loosely pinched (2) or pivotally supported (3). In this case, the compressive forces P and Q are evenly distributed 


along the edges X =+a/2 and Y=+b/2 respectively. The cases are also considered when two parallel edges of the 
plate Y =+b/2 are free from loads, and the other two are freely pinched (4) or pivotally supported (5). In these cases, 
the compressive load is applied only to the edges X =+a/2. 


Y* x 
Let us assume that F=®-—P 5 Q a Then, for the function ® , homogeneous boundary conditions are valid 


®,, =D, =0 at [X|=a/2, O,, =®,, =0 at /¥|=b/2. These conditions can be replaced by equivalent boundary 
conditions ®=®, =0 at |X|=a/2; ®=@,=0 at |Y|=b/2. Let us move on to dimensionless variables by 


formulas: 


KD, poling 


. _qD, . 
Ae ae 


b 


X=ax; Y=by; 6 » . @(X,Y)=Df (xy); W(X,Y)=w(xy)h; K= 
a 


Then system of equations (1) and boundary conditions (2), (3) can be rewritten as 


Peshkhoev I. M., et al. Critical Loads of Uniformly Compressed Orthotropic Rectangular Plate on an Elastic Base 


ees +2¢,8°O?0?w + c,0'w + 8‘kw = 8? [w,f]—p3°d2w —qd?w, é 
S*d%f +2c,02025"f + c,04f = -ac,8° [w, w]. 
[w.d,w.f,d,f]},_,, =05 [w, Aw fd, |, = (7) 
[w.dw.f.0f |, =05 [ w.2w. £0, | Lag (8) 
[w.d,w.f,d,f],_,, =95 [O,w +v,0,w,d,w +(2—v,)0,0,w, £,8,f ] =O, (9) 
[w.aw.f.aF |, =0; [@w+v,dw, dw +(2-v, )a2d,w, a = (10) 


For any values of the parameters k, p, g the nonlinear boundary value problem (6) with any of the conditions (7)— 


(10) satisfies a trivial solution (w..f,) = (0,0). 
Let us denote by E’ the set of vector functions f = (f.6,) with a finite norm, which is determined using the scalar 


product: 


(f.8) =[J(fieitfs.)dxdy; f =(f,,f,).¢=(g,.82)- (11) 
Q 
By E! we denote the space of infinitely differentiable vector functions u=(w,F), v =(w,,F) in the domain 


Q= {(x,y):|x| e1/2,, 


domain Q . We define the norm in E' using the scalar product 


(u,v). = | ie Sa (12) 
ivjs4 2 


y| <1/2} , which satisfy one of the boundary conditions (5) or (6) on the boundary OQ of the 


Ox'dy! , Ox'dy! 
The boundary value problem (6) with one of the boundary conditions (7)—(10) can be considered as a nonlinear 


operator equation 


M,u=Ih, u=(w,f)eE', (13) 
where 
5°O¢w + 2¢,8°0.0,w +¢,0)w + kw + p60. w + g0,,w a 3°[w,f] as 
— 5 U= : 
*  8tatf +. 2¢,0°628"F +0,0°F —ac,5°[w,w] 


Here M, — linear, and II — nonlinear operators, since 
[w.f] — 0. w0,f +0,w0;f —20,0,w0,0,f . 

From the results of the works of I. I. Vorovich and N. F. Morozov, it follows [9] that the operators M, and II act 
from the space E! to E’. Let k = 0, p = 0, q = 0. Then equation (9) has a unique trivial solution u, = (w.,f.) = (0,0). 
With small changes in the values of the parameters p, q and k the trivial solution remains the only one. If, for any 
small values of a number ¢>0 and the fixed value k there are parameter values p, q, that satisfy the inequality 
[P—Po|+|a—q|<« and for which there are at least two solutions to equation (9), then the point (py, qy,) is called the 


branching point (or bifurcation point) [9] of this equation. It is shown in [9] that in order to find the bifurcation point 


(Po. qo) , it is necessary to solve the eigenvalue boundary value problem, which is a linearized equation on a trivial 


solution 


Mu =0. (15) 
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The eigenvalues (Po. qo) are the critical values of the parameters of the compressive load (p. q) . Considering (13) 


and (14), equation (15) can be written as the following eigenvalue boundary value problem: to find values of the 
parameter p, such that for fixed values of the base parameter (modulus of subgrade reaction) k and the load parameter q 


along the y axis, the following problem has a nontrivial solution: 


3°Olw + 2c,8°0.0,,w + c,0,w + kw + q0,w =—ps’d-w , (16) 
[w.d,w]}, y= 9; [w.3,w =0, (17) 

[ w.a2w | pein = 83 [ w.2w | ano (18) 

[w.2,W],. = 0% [Owtvd.w,d,w+(2—-v,)OV,w, f,8,f ]) =0, (19) 
[ w. aw | wa 73 [ ow +v,d,w, dw +(2-v, )d20,w, ,8,f ] =0. (20) 


In this case, the critical value of the loss of stability of the compressed plate corresponds to the minimum eigenvalue 
of the parameter p, and the eigenfunction corresponding to this value shows the form of equilibrium after the loss of 


stability. For the case of the hinge support of the edges (18), exact solutions can be written for equation (16): 


(5‘m* +2c,m’n* +o,n*) x! —qn’n? +k ine 1 
: > Wan =Ssin7m ee sin ys |: (21) 


Pin nm 
where m and n are arbitrary natural numbers. 
From the analysis of formulas (21), it follows that with different ratios of the plate edge lengths 6, the physical 


parameters of the orthotropy of the plate material c,, c,, the parameter g of the compressive load along the y axis, and 


the modulus of subgrade reaction k (characterizing the reaction of the base), one or more half-waves in the form of 


equilibrium (21) can correspond to the critical load value p, = MIN Pian In this case, the so-called buckling by two 
proper forms is possible, when, for example, p, = MIN Py =P; =P2,- Such cases are listed below. To solve the 


problem (16) with boundary conditions of free pinching of the edges (17), a variational method is used in combination 
with a finite-difference method, which was justified in [9, 10] and is easily transferred to the case of an orthotropic 
plate. At the same time, the computer program developed for this task was tested on the solution of the 
problem (16), (18). In addition, a computer program has been developed for solving the variational method in 
combination with the finite difference method of the problem (16) with each of the boundary conditions (19), (20). For 
these cases, test examples of the corresponding boundary value problems have been constructed. 

In [7], in the problem of asymptotic integration of the equation of small oscillations of a long (b << a) rectangular 


orthotropic plate, an expansion of the solution by small parameter degrees 6=b/a is constructed and it is shown that 
the main expansion is determined from the equation of vibrations of a beam with an elastic modulus E=E,. Repeating 


the reasoning of works [7] and [10], we will look for a solution to problem (16) with boundary conditions (19) or (20) in 


the form of series: 


p=) p,8", w= ow," . (22) 
i=0 i=0 
It can be shown that the main term of the asymptotic expansion (22) of the solution to problem (16) with boundary 


conditions (19) or (20) is determined from the problem in dimensionless variables about the critical loads of a 


compressed elastic beam on an elastic base with freely pinched or movably hinged edges: 


(1-v,) wi'+ kw, +powy =0, (23) 
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WwW, =Wo =0,~ (24) 

WwW, =Wo =0. (25) 

Here the prime symbol indicates the ordinary derivative of a function depending on the variable x. It is possible to 

construct an algorithm for determining all coefficients of the terms of the expansions (22), but the authors will not do 
this in this paper. 

When moving to dimensional variables from (23)—(25), we obtain the problem of critical loads of a compressed 


beam with an elastic modulus E=E,. The eigenvalue problem (23) with boundary conditions (25) has an exact solution 


Pp, =2'm + 


a ; w,, (x) =sin(zmx). (26) 

In the case of boundary conditions (24), problem (23) is solved by the variational-difference method. 

Research Results. Numerical calculations of the critical values of the load parameter p compressing along the 
direction of the x axis with fixed values of the modulus of subgrade reaction k and fixed values of the load parameter q 
compressing along the direction of the y axis and the specified values of the orthotropy parameters of the plate material 
are carried out. The plate equilibrium forms corresponding to critical loads are constructed. 

Example 1. Let us consider an example of a square orthotropic plate from work [16] "plywood of the first type" — 


three-ply: if you direct the x axis along the fibers of the face (upper and lower of the three layers), then Young's 


modules along the fibers and across, respectively, are equal to E, =1.4-10°kg/cm* and E, = 10" kg/cm’; Poisson 


coefficients — v, =0.46 and v, =0.46/12; shear modulus for the main directions of elasticity G=1.2-10°kg/cm’. 


b D, £E 1 

Then, for the problem in dimensionless variables (16), we obtain 6 1s 2 2 0.083 ; 
a D, £E, 12 

c, =— =0.207 taking into account the above formulas for the transition to dimensionless variables. Let the edges of 


1 


the plate be pivotally supported. Then, from (21) for k = 0, q = 8.193, it follows p, = MIN Py =P, =P)» = 6.580. 


Two forms of equilibrium correspond to this critical value. Table 1 shows the values of the critical load p’ loss of 
stability of an orthotropic plate with pinched edges compressed along the x axis for different values of the modulus of 
subgrade reaction k and different values of the compressive load along the y axis. Table | shows that with an increase in 
the value of the modulus of subgrade reaction, the critical load of the loss of stability of the compressed orthotropic 
plate also increases. 

Table | 


Dependence of the critical load p’ on the base parameter k at different g 


z 
k q=0 q= 10 q=20 
0 47.455 38.011 5.298 
1,000 112.142 105.841 88.1469 
| 2,000 | 146.99 | 144.252 | 122.039 | 


3,000 164.809 163.198 141.578 
4,000 177.741 176.455 175.106 
5,000 189.002 187.843 186.65 
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Example 2. Let us consider the same orthotropic plate from Example | and direct the x axis across the fibers of the 


face. Then Young's modules along the x and y axes are E, =—-10%kg fom’ E, =1.4-10°kg/cm” respectively; 


Poisson's coefficients — v,=0.46/12 and v,=0.46; shear modulus for the main directions of elasticity 
5 2 er : : ; D,_E, 

G=1.2-10°kg/cm’ . Then, for the problem in dimensionless variables (12), we obtain 6 1; c, D E 12; 
1 1 


D 
C3 = a. = 2.481. Let the boundary conditions (13) be fulfilled — the edges of the plate are freely pinched. As a result 


1 
of numerical calculations using a computer program for solving the problem (12), (13) for this orthotropic plate, 
dependencies are constructed between the load g on the y axis and the critical value p, of the load p on the x axis at 
fixed values of the remaining parameters of this problem. Table 2 shows the values (q, p *) atk=0. 
Table 2 


Dependence between g and p* atk=0 


No. 1 2 3 4 ) 6 
q 0 100 200 300 400 500 
p* 254.388 | 241.867 | 227.789 | 208.778 | 166.462 77.105 


It should be noted that all six points (q, p*) correspond to the form of equilibrium, gradually transforming from the 
form shown in Fig. | a to the form in Fig. 1 c. Figure 1 illustrates the first, third and sixth eigenfunctions corresponding 


to the critical points (q,p*) of Table 1. These eigenfunctions show the forms of equilibrium of the plate after loss of 


stability (forms of loss of stability). 


Fig. 1. Forms of plate stability loss: a — at q = 0, p = 254.388; 
b — at q = 200, p = 227.789; c — at g = 400, p = 166.462 


As a result of numerical calculations using a computer program for solving the boundary value problem (16), (17) 


for this orthotropic plate, parameter values are found k= 200, q=0 for which two eigenfunctions w, and w, 
correspond to the eigenvalue p* = 255.6. This case corresponds to the loss of stability of the plate in two proper forms. 


Figure 2 shows graphs of the forms of stability loss corresponding to this case. 
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Fig. 2. Two forms of stability loss: a— w,; b— w, 


Example 3. Let us consider an orthotropic plate with parameters E, =1.4-10°kg/cm’, E, = 10" kg/cm’; 


Poisson coefficients — v,=0.46 and v,=0.46/12; shear modulus for the main directions of elasticity 


G =1.2-10°kg/cm’. Let the length of the plate be a=5, the width — b=1, the edges x = + — are pinched, and 


1 
the edges y= “2 are free, the parameter — 6 = 5 (the plate can be considered long), modulus of subgrade reaction — 


k =350. Numerical calculation of the solution of the problem (16) with boundary conditions (17) gives the value of the 


critical load p =37.42, the graph of the corresponding eigenfunction is shown in Fig. 3a. Numerical calculation of the 
solution to the problem (23), (24) at v, =0.46 and k=350 for a beam of length a=5 gives p, =39.01, while the 
eigenfunction has the form shown in Fig. 3b. As you can see, the critical load of the plate stability loss coincides with 


the critical load of the beam with an error |p—p,|=|37.42—39.01|~ 0.04. The eigenfunctions of problems (16), (17) 


and (23), (24) have the same number of half-waves. 


b) 


Fig. 3. Forms of stability loss: a— w (x y) , compressed orthotropic long plate; b — wy (x) , compressed beam 


Discussion and Conclusion. The problem of critical loads of an orthotropic plate compressed in two directions 
lying on an elastic base is investigated. The problem is reduced to solving a parametric linear boundary value 
eigenvalue problem. In the case of boundary conditions of a movable hinge support, exact formulas of eigenvalues and 
eigenfunctions are given, and in the case of free pinching of the edges, a variational method is applied in combination 


Mechanics 


221 


http://vestnik-donstu.ru 


222 


Advanced Engineering Research 2022. V. 22, no. 3. P. 214-223. ISSN 2687-1653 


with a finite difference method and a computer program for solving the problem is constructed. The following 
conclusions are obtained: 

1. When the compressive load component gq increases along the y axis, the critical value p* of the load compressing 
the plate along the x axis decreases. 

2. If an orthotropic plate is compressed by a load along a direction that corresponds to greater bending stiffness, then 
the critical value of the loss of stability is greater than the critical value of the compressive load acting along the 
direction of lesser bending stiffness. 

3. In the presence of an elastic base, a greater value of the modulus of subgrade reaction k > 0 corresponds to a 
greater value of the critical value of the compressive load. The presence of an elastic base increases the bearing capacity 
of the compressed plate. 

4. In the case of a long orthotropic plate on an elastic base, the critical load of a compressed plate can be calculated 
approximately as the critical load of a compressed beam on an elastic base with an elastic modulus equal to the elastic 
modulus of the plate in the longitudinal direction. 
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Abstract 

Introduction. One of the commonly used methods for assessing the dynamic characteristics of a material is the Taylor 
test, which establishes the relationship between the dynamic yield strength of a cylindrical sample material and its 
length after hitting a non-deformable barrier. The purpose of this work was to study the microhardness and determine 
the dynamic yield strength of copper samples for various impact velocities in the Taylor test. 

Materials and Methods. Experiments were carried out with cylindrical copper (M1) samples. The throwing conditions 
were selected on the ballistic stand, which provided the speed of the sample in the range of 150-450 m/s at the exit from 
the barrel. After the impact, the microhardness of the samples in the section plane was measured. The calculation of the 
dynamic yield strength was carried out according to the classical Taylor formula. 

Results. Experimental data are presented for cylindrical copper samples upon impact on a rigid wall with velocities in 
the range of 162-416 m/s, including configurations and sizes of images before and after impact. Microhardness 
distributions in the axial section of the samples were obtained. For each sample, the dependences of the averaged values 
of microhardness were constructed, which made it possible to identify four areas of deformation of the samples (the 
area of elastic deformations, plastic deformations, intense plastic deformations, the area of the material undergoing 
destruction) and determine their sizes. The dynamic yield strength of copper in the studied range of impact velocities 
was calculated. 

Discussion and Conclusions. The values of microhardness in the entire considered region and for all studied impact 
velocities exceeded the initial value. There was a significant increase in the value of the dynamic yield strength 
compared to its static value. The correlation of the maximum averaged values of microhardness and dynamic yield 


strength, which grew with increasing impact velocity, was identified. 
Keywords: Taylor test, copper cylinder, high-speed impact, microhardness, deformation, dynamic yield strength. 
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Introduction. The development in the scientific and technical fields related to the dynamic loading of bodies 
depends largely on the creation of new materials with specified properties. This was the impetus for active 
experimental, analytical and numerical studies related to the dynamic loading of deformable solids [1—4]. 

One of the commonly used methods for assessing the dynamic characteristics of a material is the Taylor test 
(method, task). The Taylor method establishes the relationship of the dynamic yield strength of the material of a 
cylindrical sample and its length after impact on a non-deformable barrier (rigid wall). This approach is often used to 
determine the dynamic yield strength of new materials [5—8], as well as to choose the determining ratios and select 
constants under the numerical modeling [9-14]. 

Light-gas guns (LGG) are used to accelerate bodies with a given shape. These installations make it possible to 
obtain a throwing speed of up to 7—9 km/s, in some experiments — up to 11 km/s, thanks to which they have been 
widely used in gas dynamics, ballistics, materials science, etc. The Research Institute of Applied Mathematics and 
Mechanics, Tomsk State University, has developed a single-stage light-gas gun [15] in which the sample is accelerated 
by compressed gas (helium) supplied from a balloon. It was used to conduct experiments by the Taylor method 
presented in this paper. 

This paper is aimed at studying the microhardness and determining the dynamic yield strength of copper samples for 
various impact velocities in the Taylor test. 

Materials and Methods. The experiment was carried out with cylindrical copper (M1) samples with a length of 
34.5 mm, a diameter of 7.8 mm, and a weight of about 15 g. The composition of the sample material is indicated in 
Table 1. 

Table | 
Composition of copper (M1) 


Cu Ag Fe Ni S AS Pb O Sb Bi Sn 


>99.9 | <0.003 | <0.005 | <0.002 | <0.004 |<0.002 | <0.005 | <0.05 |<0.002 | <0.001 |< 0.002 


The throwing conditions were selected on the ballistic stand, which provided the speed of the sample in the range of 
150-450 m/s at the exit from the barrel. The selection of conditions was required, since samples of significantly smaller 
mass were usually used for such experiments, which were thrown at a much higher initial velocity. By adapting the 
initial conditions and equipment, it was possible to provide stable throwing of the copper cylinder at a given speed. 
After the experiment, the samples were cut into two equal parts along the axis of symmetry using DK7732 CNC 
machine for electroerosion cutting. 

Microhardness was measured in the section plane along the axial line of the samples on a PTM-—3 hardness tester 
according to GOST 9650-76 by indentation of diamond tips. The measurement error of this device was 2 %. The 
microhardness values were also calculated according to GOST 9650-76. 


The calculation of the dynamic yield strength o was carried out according to the classical Taylor formula: 
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where p — density of the material; » — impact velocity; Ly) — initial length; L — final length after impact; L, — length 
of the elastic part of the sample. In this formula, the key indicators were the length of the elastic part of the sample and 
the final length of the cylinder after impact. 

Research Results. Microhardness. Figure 1 shows cross sections of copper samples after impact at different initial 
velocities. For all samples, there was a region of elastic deformation with a finite diameter equal to the initial one, 
including for a speed of 416 m/s. The elastic region turned into a plastic one, which was accompanied by deformation, 
including in the radial direction, and, accordingly, an increase in the final diameter. Closer to the contact surface, a zone 
of intense plastic deformations that passed into the zone of destruction of the cylinder material was observed. There was 
also a slight asymmetry of the deformation of the samples due to the characteristic properties of throwing, which might 
create some difficulties for direct comparison of the results of experiments and numerical modeling when using the axis 


of symmetry or the plane of symmetry in the formulation of the problem. 


Ba 


c) d) 


Fig. 1. Cross sections of copper samples after impact at different initial velocities: 
a — 162 m/s; b — 225 m/s; c — 316 m/s; d— 416m 


a) 
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Fig. 2. Distribution of values along the centerline of samples for different impact velocities 
(1— 162 m/s, 2— 225 m/s, 3— 316 m/s, 4— 416 m/s): a — microhardness; 
b — calculated averaged distribution of microhardness 


In the initial state, the measured average microhardness was 1,150 + 100 MPa. Such a spread of values, apparently, 


occurred due to the structural features of the material and the presence of inhomogeneities in the structure. Figure 2 a 
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shows the distribution of microhardness values along the centerline of the sample for different impact velocities, and 
Figure 2 b presents the results of averaging the obtained values. 

Through comparing the average values of microhardness to the value in the nondeformable sample (1,150 MPa), we 
obtained that microhardness in the entire region under consideration exceeded the initial value. Several deformation 
zones can be distinguished. This is the rear region of the sample, where there is a slight increase in microhardness to 
values of 1,250—1,350 MPa, apparently associated with the effect of gas pressure during throwing. In the middle part, 
microhardness approaches the indicators of undeformed samples. Further, closer to the contact boundary, 
microhardness starts to grow to values of 1,400-—1,600 MPa with an exit to the inflection point, after which there is a 
sharp increase in microhardness to 1,800—2,700 MPa. 

The maximum microhardness is observed at an impact velocity equal to 316 m/s, while at higher speeds, a sharp 
decrease in the microhardness value is observed in Figure 2. The drop in microhardness in the sample for the impact 
velocity of 416 m/s is due to its destruction in the impact area and the loss of part of the cylinder material in the form of 
fragments. 

Let us take a closer look at a sample with an impact velocity of 316 m/s. For this sample, two series of 
measurements were made along two lines, which are notionally named C and Cl. The location of these lines was 
chosen from the following considerations. Line C was located along the axis of symmetry, and line Cl passed through 
the middle of the radius of the sample section. During the microhardness measurements, over 100 measurements were 


made along each line. This array of microhardness values was averaged, the result of averaging is shown in Figure 3. 
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Fig. 3. Averaged distribution of microhardness of a copper sample at an impact 
velocity of 316 m/s along the lines: a— C; b— Cl 


Figure 3 shows that both curves have a nonlinear character, and the distribution of microhardness in both cases is 
qualitatively similar and quantitatively close. Such data make it possible to identify areas of deformation of samples and 
determine their dimensions. 

Areas of deformation of cylindrical samples. It is proposed to consider the scheme of deformation of a cylindrical 
sample, shown in Figure 4 and including four areas, whose size is determined based on the analysis of the distribution 
of microhardness. 

Zone | corresponds to the area of elastic deformations; 2 — plastic; 3 — intense plastic deformations; 4 — the area 
of the material undergoing destruction. Table | shows the sizes of these sample areas depending on the impact velocity, 


where v — impact velocity; L — final length after impact; L, — length of the elastic deformation zone; L, — length of 
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the plastic deformation zone; L;,;— length of the zone of intense plastic deformations; L;— length of the fracture zone; 


D, — diameter of the rear end of the cylinder; D, — diameter of the contact boundary. 


Di 


D2 


Zone | 


Zone 2 


Zone 3 
Zone 4 


Fig. 4. Scheme of sample deformation after collision with a rigid wall 


Geometric dimensions of the deformation areas of the samples 


v, m/s Lmm L.,.mm L,,mm. Lip, 1AM L;mm D,,mm Dz,mm 
0 34.5 0 0 0 0 78 78 
162 26.1 12.34 10.63 3.13 0 7.8 12.7 
225 22. 9.78 9.51 3.21 0.93 7.8 15.8 
316 16.1 5.87 6.83 3.4 1.14 79 21.4 

416 9.3 3.5 1.38 4.42 2.01 79 31.28 


Table 2 


The data given in Table 2 show that the sample had no material destruction zone at an impact velocity of 162 m/s. 


The resulting sample is an example of the classical Taylor test, which can later be used to create an adequate numerical 


model of the impact of a cylindrical sample on a rigid barrier, and the selection of constants of the material models. 


Samples at speeds of 225 m/s and 316 m/s exhibited all four zones, but compared to the sample at a speed of impact of 


416 m/s, the destruction zone was small. A cylinder with an impact velocity of 416 m/s had practically no plastic 


deformation zone after the test. At a given impact velocity, the elastic deformation zone quickly turned into a zone of 


intense plastic deformations combined with a fracture zone. 


Dynamic yield strength. The results of the calculation of the dynamic yield strength o are presented in Figure 5 a. 


The speed value of 416 m/s was not taken into account, since in this case, the destruction of the cylinder was 


significant, which did not allow applying this calculation method. The static yield strength of M1 copper was 0.1 HPa. 


There was a significant increase in the value of the dynamic yield strength compared to its static value; there was also 


an increase in the dynamic yield strength with an increase in the impact velocity. 


Pakhnutova N. V., et al. Microhardness and Dynamic Yield Strength of Copper Samples upon Impact on a Rigid Wall 


o, hPa Hm, MPa 
0.44 - 
ee 2,500 4  a—-_£ 
0.3 5 ° 2,000 4 
0.2 5 1,5007 
‘i 
0.14 ; : ; = 1,000 4 1 1 T 
0 100 200 300 vu, m/s 0 100 200 300 ~—sv, m/s 
a) b) 


Fig. 5. Dependences on the impact velocity: a — dynamic yield strength; 
b — maximum microhardness 


Figure 5 b shows the dependence of the maximum microhardness of the samples on the impact velocity. The graphs 
in Figure 5 a and 5 b are qualitatively similar, which implies a correlation of the dynamic yield strength and maximum 
microhardness, and the possibility of their mutual recalculation. 

Discussion and Conclusions. The results of Taylor test experiments for cylindrical copper samples in the range of 
impact velocities of 162-416 m/s are presented. The data obtained has shown that after the impact of the cylinder on the 
rigid wall, the microhardness in the entire sample exceeded the initial value of 1,150 MPa, and in the impact area there 
was a significant increase in microhardness up to 1,800—2,700 MPa. The separation of the deformed cylinder into four 
areas was proposed: elastic deformation, plastic, intensive plastic deformation and destruction. Estimates of the sizes of 
these areas for the studied impact velocities were given. According to the classical Taylor formula, the dynamic yield 
strength was calculated, which significantly exceeded the static yield strength and grew with increasing impact velocity. 
The dependences of the dynamic yield strength and the maximum averaged value of microhardness on the impact 
velocity of the sample on a nondeformable barrier were shown. The correlation of the maximum averaged values of 
microhardness and dynamic yield strength, which increased with growing impact velocity, was specified. The presented 
data can be of value for assessing the adequacy of the physicomathematical model used for numerical calculation of 


problems of high-speed deformation of metals and alloys. 
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Abstract 
Introduction. In sheet and hull structures operating under pressure, destruction, as a rule, is localized along the transi- 


tion line from the base metal to the weld metal. Methods of increasing the durability of butt-welded joints, which are 
aimed at reducing stress concentration and creating favorable residual compression stresses, are described. 
Materials and Methods. The tests were carried out on an installation for biaxial bending, which created a biaxial stress 
field. Factory-made coupons and samples with an additionally processed transition zone from the weld metal to the base 
metal were tested. The effectiveness of further processing is shown by the following methods: 

— abrading; 

— grit hardening; 

— abrading with grit hardening; 

— melting of the fusion line in argon without filler wire; 

— melting of the fusion line in argon with filler wire EP-410U; 

— melting of the fusion line without filler wire with plastic deformation between narrow rollers. 
Results. The origin, development of destruction, and its features were analyzed using different methods of further pro- 
cessing of joint welds. Confidence spans (95 %) of the origin and development of failures for joint welds and base metal 
were calculated. The efficiency of the proposed methods for further processing was evaluated. 
Discussion and Conclusions. An analysis of the effectiveness of methods for increasing the durability of butt-welded 
joints has shown that the creation of a smooth transition from the weld metal to the base metal reduces significantly the 
stress concentration. This provides increasing the number of cycles before the onset of destruction and the survivability 
of compounds. Due to compressive stresses in the near-weld area, it is possible to increase the durability of joint welds. 
The most effective methods of further processing of welds combine the reduction of stress concentration and the crea- 
tion of residual compression stresses. The high-tech solution is remelting the transition zone in an argon medium with 


an additional EP-410U filler wire. 


Keywords: durability increase, cyclic loads, butt joints, stress concentration, residual stresses, origin and development 


of fracture. 
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Introduction. Welded hull and shell structures operating under pressure are widely used in shipbuilding, chemical, 
power engineering and other industries. Reducing the metal consumption of such structures while increasing their effi- 
ciency is a critical task [1-4]. To reduce the weight of structures, it is required to use stronger materials under stresses 
close to the yield point. In this case, special requirements are put forward for the implementation of welded structures to 
ensure high structural strength. In other words, it is required to develop measures that increase the durability of struc- 
tures to values close to the indicators of the base metal. 

The work is aimed at studying the possibilities of increasing the durability of welded joints operating under cyclic 
loads in a corrosive environment (3 % NaCl solution), with the help of additional processing of joints. 

Materials and Methods. Butt-welded joints were made of chromium-nickel-molybdenum steel. Preliminary heat 
treatment of sheet material provided yield limits of 900 MPa; 1,100 MPa; 1,150 MPa. The tests [5, 6] were carried out 
with a biaxial stress field and simultaneous action of repeated static loads in a corrosive environment (3 % NaCl solu- 
tion). Manual multi-pass welding was performed with low-alloy electrodes 48H 11, 48H13, and austenitic electrodes 
EA 981/15. The failure nucleated from the stretched fibers, so the stress state was studied on the stretched surface of the 
sample. The stresses were determined by calculation. They were measured with a lever Huggenberger tensometer and 
tensoresistors with a base of 5 mm at a distance of 10 mm from the weld. If the measured stress values differed from the 
calculated values by more than 5 %, the pressure under the sample was corrected. 

The failure of welded joints was localized along the transition line from the base metal to the weld metal, as shown 


in Figure 1. 


Fig. 1. Fracture pattern of the butt-welded joint (photo of the authors) 


The main reasons for the decrease in the working capacity of welded joints compared to sheet metal [7—11]: 

— geometric stress concentration; 

— residual welding stresses; 

— deterioration of the structure and properties of the base metal in the near-weld area under the influence of the 
thermal welding cycle. 


The influence of each of these factors can be partially or completely neutralized. 
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At the Department of “Machines and Automation of Welding Production”, DSTU, butt-welded joints of high 
strength steels have been tested by various post-treatment methods for many years (Table 1). The experiments corre- 
sponded to the loading conditions of real structures. 

Table | 


Techniques of additional treatments of butt-welded joints in the transition 
zone from the joint to the base metal 


No Treatment 
1 | Cleaning with an abrasive wheel with a grain size of 80. Fillet radius R = 30-40 mm 
Grit hardening with AD-1 shot blaster with DSL—1.5cast steel shot. Pressure is 5 atm, head travel speed is 
75 mm/min. Processing zone at the transition point from joint to metal is 15-30 mm 
Cleaning the weld metal — base metal transition zone with an abrasive wheel by the first method, and grit harden- 


ing — by the second method 
Washing in argon medium without filler wire. Diameter of the tungsten electrode — 3 mm, current 
4 |1=120A, are voltage U = 10-12 B, welding speed is 8 m/h, transverse vibration frequency — 60 min’, vi- 


bration amplitude — 6 mm 


5 | Washing with EP—410U filler wire by the fourth method. Filler wire diameter is 1.6 mm 


Washing without filler wire by the fourth method. Plastic deformation between the beads. Bead diameter — 
120 mm, width — 20 mm. Bead pressure — 18,000 kgf, rolling speed — 1.4 m/min. 


Physical and metallurgical processes occurring under welding cannot be modeled in full. Therefore, the tests were 
carried out on full-scale butt-welded joints with full preservation of the factory welding technology and geometric pa- 
rameters of the joints. Low-cycle fatigue of welded joints was studied on samples in the form of disks with a diameter 
of 550 mm and a thickness of 30 mm on UDI-550 installation [12]. They were pivotally fixed along the contour and 
loaded with hydrostatic oil pressure. Under the action of hydrostatic pressure, the disk was axially bent. Compression 
stresses appeared on the inner surface, tensile stresses appeared on the outer surface, and a corrosive medium, a 3 % 
aqueous solution of sodium chloride, affected it. The samples were tested under repeated static loading with a frequency 
of 10 cycles/min. 

The maximum stresses occurred in the central part of the samples. On a large surface, the statistical probability of 
the occurrence and development of failure is higher, which in general brings the test conditions closer to the real work- 
ing conditions of the loaded hull structures. 

Research Results. Table 2 shows the test results of samples after additional processing of welded joints under cy- 


clic loading. 


Table 2 
Influence of additional methods of processing butt-welded joints 
on performance characteristics 
. Additio- Number of cycles 
Yield : a Max cycle before 
Welding Type of joint nal . . 
No. | strength, ; . voltage, appearance | loss of Breaking point 
materials processing process- . 
MPa MPa . of cracks, | tightness, 
. Nt Np 
1 No 2,600 20,560 
2 EA lst method: abrad- 3,890 24,401 a ; 
—— 1,1 7 ———————_ ———————COCWrSTT tion line* 
3] 110° | oginis ing oy . 10,840 | 27,237 (nrecle arena 
mn “ 9115 | 28,947 
5 baa seetiodi oat N 2,000 7,430 Transition line 
: (o) 
6 | 1,00 | agmi3 | “Seon 860 2,100 | 11,500 
hardening er 
7 Yes 1,410 27,349 Transition line 
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. Additio- Number of cycles 
Yield : 8 Max cycle before 
Welding Type of joint nal : ; 
No. | strength, sentra ene voltage, rocess. | aPPearance loss of Breaking point 
MPa P g MPa P ” of cracks, | tightness, 
a Nt Np 
8 2,000 14,051 
9 Y 10,670 32,310 
1,100 48H13 | 3rd method: abrad- 700 . Transition line 
10 P . 12,430 37,540 . 
ing and grit harden- with access to the 
11 1.100 EA a” 720 Yes 10,870 32,840 ee 
12 : 981/15 . 9,300 34,460 
13 3,800 24,930 
No 
si 48H13 | 4th method: h aka mae T. ition li 
5 aa te) vas - Yes 6,240 33,287 —— ion line 
1,100 ing in argon without 700 with access to the 
16 . 6,450 29,714 
filler wire base metal 
17 EA Yes 8,320 25,400 
18 981/15 7,000 24,250 
19 4,300 29,074 i . 
0 No 4,200 33,079 Transition line 
. 900 | 48H11 | a. i = 605 | a 
2 ee ae sas 7 29,800 | 58,384 | aes oe 
t 
3 wi wire es 30,910 91,300 and perpendicular to 
the weld 
24 24,100 86,000 
= No aa ae Transition line 
1,1 48H11 th method 72: : : 
oy) eee eiumene a : 7,627 | 39,457 Transition line 
28 - 4,890 | 33,405 and base metal 
= No a oa Transition line 
1,150 48H13 5th method 760 : : — - 
31 Y 8,500 36,400 Transition line 
32 ra 9,886 | 34,636 and base metal 
33 3,270 21,980 
oa No 1300 oo Transition line 
35 1,200 19,079 
36 3,160 26,880 
48H13 605 
37 20,860 48,210 B Sand 
- Washing of fusi Yes ea Eas hee. 
39 | 900 ee ee 17,300 | 93,552 
line weld 
40 ‘ 24,475 99,910 
41 22 
a 6th method: weld No aa — Transition line 
9A981/1| rolling between : ; 
43 5 narrow beads 720 31,950 96,875 Base metal and 
Yes cracks across the 
44 23,450 78,543 weld 
45 No 3,285 23,245 Transition line 
46 | 1150 | 48H13 760 22,900" |), -87,280%.|  -Bavemeial and 
Yes cracks across the 
47 21,168 55,039 weld 
*from weld to base metal. 


Additional processing of welded joints by the Ist and 4th method (Table 1) reduced the stress concentration through 


increasing the interface radius of the weld metal and the base metal. 
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The 2nd method created small compressive stresses in the weld area, but practically did not change their concentra- 
tion. An alternative to the proposed method is rolling the transition zone, presented in [13]. 

The third, fifth and sixth methods, in addition to reducing the stress concentration, made it possible to obtain favora- 
ble residual compression stresses in the transition zone from the weld metal to the base metal. 

Figure 2 a shows the confidence intervals (95 %) of the failure initiation, and Figure 2 b — intervals before failure 
for welded joints (dotted lines) made according to the factory technology, and for the base metal (solid lines). The in- 


tervals were plotted according to the data from Table 2. 
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Fig. 2. Comparison of the operability of the base metal and butt-welded joints under repeated static loading: 
a — number of cycles before appearance of cracks; b — number of cycles before failure 


The spread of values of resistance to damage and durability of butt-welded joints (Fig. 2a) and the base metal 
(Fig. 2 b) depending on the maximum stresses at zero pulsating loading cycle are presented. 

It can be seen that the resistance to the failure nucleation N, and the durability of butt-welded joints before failure N, 
(in this case, before the loss of tightness) is much less than that of the base metal. 

To determine the compressive stresses generated in the joints as a result of processing, residual stresses were meas- 
ured. This method is described in [14, 15]. 

Table 3 presents the results of measurements of residual stresses in the direction perpendicular to the axis of the 
weld on the surface of the welded joints in the near-weld zone. 

Table 3 
Residual compression stresses in the interface zone of weld metal and base metal, 
depending on the additional processing methods 


Method of additional processing of the weld | Compressive stresses, MPa 
2nd and 3rd 60-80 
Sth 240-320 
6th 700-800 


According to Table 2 and Figure 3, it is possible to judge the effectiveness of various methods of additional pro- 
cessing of welded joints operating under cyclic loading. 


io} 


a | 
MPa }, 
_— Base 
900K, KK 
goo} oA 
700 My [OO] 
600 Ts y 
ss ZB 
Joint Ts ‘ 
500 | welds aN 
* 3 
400 : N,x10 
2 4 6 8 


Lyudmirsky Y. G., et al. Increasing the Durability of Butt-Welded Joints Operating under Cyclic Loads in a Biaxial Stress Field 


1st method 
2nd method 
3rd method 
4th method 
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6th method 


xX HePOO xX 


b) 
Fig. 3. Efficiency of methods of additional processing of the weld metal — base metal interface zone: 


a — number of cycles before cracks appear; b — number of cycles before failure 


Figure 3 shows the results of testing samples of butt-welded joints processed in different ways. Figure 3 a shows the 
confidence areas of the spread of values before the initiation of the destruction of the base metal and butt-welded joints 
without additional processing. Figure 3 b shows the same areas before the destruction of the samples (before loss of 
tightness). 

Dressing the weld metal — base metal transition zone increased the resistance to the fracture nucleation due to a de- 
crease in the concentration of stresses (Table 2). However, this method did not practically affect the fracture geometry. 
Cracks nucleated and developed along the interface line of the weld metal and the base metal. 

Samples had a similar fracture pattern, in which the fusion line of the weld and base metal was washed by an arc in 
argon without a filler metal (4th method). This caused a decrease in stress concentration. As a result, the resistance to 
the failure nucleation increased. At the same time, the resistance to the development of fracture did not practically 
change (Fig. 3). The 4th method is more technologically advanced compared to the Ist one, since it does not require 
additional equipment other than welding. 

Grit-hardening (2nd method) did not practically affect the damage resistance, determined by the number of cycles 
before the appearance of a visible crack — Nr, but increased the survivability — the number of cycles that the sample 
withstood after the formation of the first crack until it lost its bearing capacity (leak). This is due to the fact that grit- 
hardening does not guarantee uniformity of the surface deformation of the metal, specifically, at the junctions with un- 
dercuts, surges, craters, and non-melting outlines of the weld. It is here that failure starts. However, in places of smooth 
coupling of the weld metal and base metal, shot blasting, which caused the greatest compression stresses, prevented the 
formation of an extended main crack, which increased the resistance to the development of destruction. 

To increase the efficiency of the grit-hardening, the 3rd method was proposed. Abrading was performed by the Ist 
method, and then shot blasting followed. In comparison to the 2nd method, the durability of the samples increased by 
about 20 % before the onset of fracture and before the loss of tightness. Comparison of the Ist and 3rd methods showed 
that the number of cycles before the fracture nucleation did not actually change, but the number of cycles before frac- 
ture (loss of tightness) increased by 20 %. 

The 5th method [16] was washing with an EP-410U filler wire with a diameter of 1.6 mm. When cooled (140 °C 
and below), the fillet beads underwent martensitic transformations. When cooled to room temperature, the volume in- 
creased by 1.5 % in total [11]. As it was shown earlier, this causes the appearance of residual compression stresses up to 
300 MPa. Fillet beads with an increased specific volume contributed to an increase in damage resistance and the devel- 
opment of failure, i.e., an increase in survivability. In such joints, the first cracks appeared either on the base metal, 
away from the weld, or simultaneously along the transition line from the weld to the base metal (Table 2). Cracks found 
on the fusion line, as a rule, developed at a low rate, failure was caused by crack coalescence in the base metal and in 
the welded joint. The typical appearance of the welded joints, additionally processed by the 5th method, after the test is 


shown in Figure 4. 
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Fig. 4. Destruction of a butt-welded joint made by the 5th method (photo of the authors) 


The considered type of treatment reduces the stress concentration along the entire length of the weld and creates fa- 
vorable residual compression stresses in the area of the fillet beads with an increased specific volume. 

The action of compressive stresses inhibits the fracture nucleation. With the development of cracks along the fillet 
bead, the action of transverse compressive stresses is no longer effective and is partly eliminated. This can explain the 
absence of the influence of beads of increased specific volume on the survivability of welded joints. If the failure nucle- 
ates not on the fillet bead, but on the base metal, the survivability of welded joints increases. 

Treatment of the transition zone from the base metal to the weld metal by the 5th method increased the resistance to 
failure nucleation and development approximately 3-4 times. At the same time, the topography of the failure changed 
markedly. The first cracks, as a rule, originated in the base metal (Fig. 4). In the presence of poorly welded craters on 
the front surface of the weld, cracks initiated in these places and developed mainly perpendicular to the weld, but even 
in this case, the resistance to the failure nucleation and development remained quite high (Table 2, 5th method). 

Surface plastic deformation of the near-weld zone [12] in butt joints by rolling with narrow beads (6th method) 
made it possible to bring the resistance to nucleation and development of failure to the level of similar characteristics of 
the base metal. This conclusion was confirmed by the fracture pattern of welded joints made by the 6th method (Fig. 5). 
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b) 
Fig. 5. Fracture nucleation and development pattern of the welded joint, 
whose interface zone was run between narrow beads (6th method): 
a — record of the sample fracture development; b — appearance of the welded joint during fracture (photo of the authors) 

It can be seen that the fracture nucleated and developed according to the basic method (Fig. 5 a, Table 2, and the 
fracture pattern of the joints, shown in Fig. 5 b). The washing of the fusion line and its subsequent running between the 
narrow beads (6th method) increased the resistance to the fracture nucleation by about eight times, the survivability of 
welded joints — by about four times. 

Discussion and Conclusions 

1. As initial samples, we considered welded joints made according to the factory technology without additional pro- 
cessing of the transition zone. In this case, during cyclic loading along the transition line from the weld metal to the 
base metal, multistage nucleation of fatigue cracks was observed. They developed rapidly and were combined into one 
mainline. Then, it developed in depth, which caused loss of tightness. The durability of welded joints turned out to be 
2-3 times less than that of the base metal. 

2. Two types of additional processing were involved: 

— abrading of the transition zone from the weld to the base metal (1st method); 

— remelting of the transition zone from the weld to the base metal with a non-melting electrode, an arc burning in an 
argon medium (4th method). 

This made it possible to increase the resistance to the fracture nucleation by almost 2 times. The survivability of the- 
se welded joints did not change much. 

3. Grit hardening (2nd method) did not practically affect the number of cycles before the fracture nucleation, but 
slightly increased the number of cycles before fracture. 

4. Preliminary abrading of the weld metal — base metal transition zone and subsequent grit hardening (3rd method) 
increased the resistance of welded joints to the fracture nucleation and development by almost 1.5 times. 

5. When the fillet bead was surfaced with a material with suitable volumetric changes (Sth method), cracks occurred 
in the transition zone and in the base metal. At the same time, the resistance to the fracture nucleation and development 


increased approximately 3—4 times. 
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6. Washing without filler wire (4th method) and plastic deformation between narrow beads (6th method) provided 
an increase in the resistance to fracture and survivability of welded joints almost to the level of similar characteristics of 


the base metal. 
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Abstract 

Introduction. With the growth in the production and transportation of gas and oil, the urgency of the problem of 
tightness of the connection of pipes of the oil and gas range increases. The most common are coupling and streamline 
threaded connections. Threads of the buttress type, conical trapezoidal, and triangular threads are used. The tightness of 
the connection depends directly on their quality. The production of pipes and couplings is influenced by many factors. 
Examples include technological heredity, rigidity of technological equipment, its setting and resetting. The required 
connection quality is provided by the quality of the elements. The assembly should take into account the possibility of 
complete interchangeability. The second recommended option is a group assembling. In all cases, control operations are 
mandatory. The paper presents the results of studies on the distribution of pipe and coupling thread quality parameters 
in the batch. 

Materials and Methods. Fragments (tubular parts) of casing pipes with trapezoidal thread and increased 
tightness (CPTT) and with trapezoidal thread of the buttress profile were studied. They were randomly selected from 
different batches. When measuring, the samples were fixed in the chuck of 11611P lathe. The measuring instrument was 
a clock type 0 accuracy class indicator with a stand. Generally accepted statistical methods were used to evaluate the 
data obtained. 

Results. The experiments have found how the tightness of the connection through the method of complete 
interchangeability is affected by the radial runout of the surfaces of the thread profile tops, the pipe chamfers, and the 
sealing belt. 

Discussion and Conclusions. Significant deviations of the values of the thread parameters from the permissible ones 
were found, which might cause a leak in the connection. Obtaining a high-quality and reliable connection is possible 
when assembling by the method of incomplete or group interchangeability. 


Keywords: OCTG connection, connection quality, quality control, buttress connection, CPTT connection, radial runout 
of buttress thread surfaces, radial runout of CPTT thread surfaces. 
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Introduction. Exaggeration of the operating conditions for pipes is associated with the complication of profiles and 
an increase in the depths of wells in new oil-and-gas fields. Oil and gas grade pipes are exposed to significant 
alternating loads at the operating temperature ranges from minus 60 °C to plus 200 °C [1-6]. The products in question 


must meet the requirements for resistance to fatigue and brittle fracture, as well as corrosion resistance in aggressive 
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environments. This determines the choice of materials for their production. The considerable length of oil and gas 
pipelines implies providing reliable connection of pipes. Domestic manufacturers offer coupling and inserted 
connections with conical trapezoidal and triangular threads [7-10]. In the coupling types, the following are 
distinguished: 

—CPTT (for casing pipes with trapezoidal thread and increased tightness); 

— with buttress trapezoidal thread. 

The connection of the first type is characterized by high tightness and resistance to tensile loads. These qualities are 


provided by the shape of the thread profile — an uneven trapezoid with a pitch of 5.08 mm and a taper of 1:16 (Fig. 1) 
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Fig. 1. Coupling connection of casing pipes of CPTT type 
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The angle of the stab flank of the pipe of 10 ° provides an easy fit into the coupling and reduces the probability of 
thread jamming. Under tensile and bending loads, the pipe thread should not disengage from the coupling thread; 
therefore, the support end of the profile is made at an angle of 3 ° [8-12]. The mating conical sealing surfaces located 
behind the thread increase the tightness. When connecting, the amount of tension is limited due to the contact of the 
products along the inner bearing faces. The sealing ring made of PTFE increases the tightness of the joint [13]. 


The scheme of the coupling connection with a trapezoidal buttress thread is shown in Figure 2. 
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Fig. 2. Coupling connection of buttress casing pipes 


A distinctive feature of this type of connection is the ability to take significant tensile loads for a long time. The 
buttress thread profile has the shape of an unequal trapezoid with a pitch of 5.08 mm and a taper of 1:16 [14]. 
The angles of inclination of the embedded and support ends are 10° and 3°, respectively. This simplifies assembly 


and reduces the probability of jamming. For tightness, a special lubricant or fluoroplastic rings are used [15]. 
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Basic requirements for common type of casing pipe connections: 

— permeability in wellbores of varying complexity and depth, including in areas of significant curvature; 

— tightness of pipe columns and high strength for all types of loads; 

— providing the passage of devices and tools into the inner space of the pipe columns during technological 
operations. 

It is required to strictly comply with the requirements of standards for the quality of threaded connections, as well as 
for the means and methods of control. The characteristics of threaded connections of casing pipes are regulated by 
GOST 632-80 and API 5B international standards. These documents prescribe to control the following parameters: 
appearance of the thread, thread geometry and sealing conical surfaces, tension on the threaded working gage. 

Note that the shapes of the mating elements of the pipe and coupling in cross sections do not always coincide. In 
these cases, in some areas, the mating surfaces will not contact, that is, the mating will remain loose. This creates 
conditions for a leak in the connection. If such a gap is not eliminated by positive drawing-up of the conical joints, then 
it is impossible to provide the effect of a ring lock, which guarantees tightness. Thus, deviations from the taper create 


prerequisites for leaks. As an example, consider the connections of CPTT and buttress pipes with an outer diameter of 


+0.45 


+0.3 
92 mm and “,,; mm. 


127 mm. On a section of 25.4 mm, the permissible deviations from the taper are respectively: 


Materials and Methods. The data obtained in the framework of the study were statistically processed using 
generally accepted methods. 

To determine the degree of tightness of the suitable pipes, the radial runout of the conical and cylindrical surfaces of 
the pipes was measured. These elements were fixed in the outside jaws of a three-lobe lathe chuck (i.e., along the pipe 
base surface on 11611P machine), subject to minimization of elastic deformations of the object from the clamping 
forces. First, the lathe jaws were brought into contact with the nozzle with a force sufficient to hang it off — this is how 
its axial position in the chuck was corrected when the spindle was turned. The surface runout was visually controlled. 
Then, in the plane of the jaws, a contact was created between the surface of the tubular part and the dial gauge mounted 
on the rack (Fig. 3). The tubular part was clamped in the chuck jaws. The gauge arrow was used to make sure that the 


clamping force did not deform the wall of the tubular part between the chuck jaws. 


Fig. 3. Scheme for measuring radial runout of the tubular parts and couplings: 


1 — product; 2 — three-lobe lathe chuck; 3 — dial gauge 


To measure the runout of the seal band and chamfer, the lathe slide with the gage holder initially remained 
stationary, and the ball attachment of the gage was installed in the middle of the band and in the middle of the chamfer. 
The gage readings were taken and recorded by turning the machine spindle every 10° (fixed in 36 sections, starting 
from zero). The gage moved in the axial direction of the part, the ball attachment maintained its position relative to the 
center of the chamfer or band. 

The runout of the threaded surface was measured along the thread ledge ribbon. To do this, the gear and the feed 
box of the machine were tuned to an inch thread with five threads per inch. After turning on the sliding clasp nut 


through spinning the spindle, all backlashes in the kinematics of the machine were selected. At that, the gage was set to 
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the middle of the ribbon of the first thread in the zero section. 


Research Results. Figure 4 shows the results of measuring the radial runout of the outer surfaces of the tubular 


parts with various connections in the polar coordinate system. 
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Fig. 4. Radial runout of the outer surfaces of threaded pipes: a — buttress; b — CPTT 
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To control the taper in the tolerance field, the deviations of the thread heads on five threads from a given profile 
(Archimedean spiral) in the radial direction were measured. For the buttress tubular part, a graph was plotted in the 
polar coordinate system based on the values of the measured radial runout of the surface of the thread heads (Fig. 5 a). 
Here, the radial scale had values from 0 to 2 mm, and the circular scale — from 0 ° to 360 °. Similarly, the radial runout 
of the thread apexes on the CPTT tubular parts was investigated (Fig. 5 b). 
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Fig. 5. Radial runout of the surface of the thread heads from the start of thread at five turns: a — buttress; b— CPTT 
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The diagram (Fig. 6 a) is based on the results of measuring the radial runout of the chamfers of tubular parts with 
buttress thread. The values of the radial runout of the band are reflected in graph 6 b, which shows the circularity 
deviation of the seal band. Comparison of the values of the circularity deviation of the band shape obtained by different 


methods shows their identity. 
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248 Fig. 6. Radial runout: a — of the tubular part chamfer with buttress thread; b — of the seal band of the CPTT connection 
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Discussion and Conclusions. The analysis of Figure 4 has shown that pipes with different types of connections are 
characterized by significant taper deviations. At the same time, the surface of the pipe with a buttress thread is relatively 
smooth. If the diagram visualizes the CPTT thread, then there are sharp outliers of a large number of points. This means 
that there is a rough microrelief on the surface of the pipe with significant random surface level differences. 

In Figure 5 a, the irregularity of the ascending spiral line with a difference of radius vectors of more than 0.4 mm is 
clearly visible. The wave-like change in the ordinate line of the profile is a consequence of the mismatch of the rotation 
axis of the tubular parts and the thread axis. This harmonic component can be compensated by special mathematical 
processing, which was not planned in this study. The analysis of Figure 5 b shows rather large irregularity of the 
ascending line of the thread head. This indicates significant deviations of the surface under consideration from a given 
profile (Archimedean spiral) for the CPTT connection. 

The analysis of Figure 6 has shown that the deviation from sphericity (difference between the largest and smallest 
radii-runout vectors) is about 0.25 mm. The shape errors of the seal band for the CPTT connection were determined by 
different methods. According to the measured values of the band diameter (in the plane of the end face), a variation in 
size up to 0.25 mm (from the zero section, through 30 °) was detected: 169.35; 169.20; 169.10; 169.15; 169.25; 169.30. 
The angular distribution of the diameter counts indicated the oval shape of the band surface in cross section. 

Thus, it can be concluded that there are significant deviations on the surfaces of pipes and couplings. We are talking 
about the quality of threads, chamfers, bands, which determines the connection reliability. These deviations cause the 
appearance of defective products in the batch. Their use does not provide the tightness and reliability of the pipe — 
coupling connection, which causes leaks during hydraulic tests. The specified connection quality is provided only by the 
method of group interchangeability, which requires preliminary separation of the assembled products into groups and 
generates additional costs. It is required to tighten the tolerances for the manufacture of threaded parts of pipes and 
couplings, which will allow for the assembly by the method of complete interchangeability. The considered deviations 
may occur due to the instability of the technological process of obtaining threads. In particular, the insufficient rigidity 


of the technological system and its untimely adjustment affects negatively. 
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Abstract 
Introduction. Pressure-operated thick-walled hull structures are the most common type of high-duty welded structures. 


When these structures are loaded with internal pressure, a complex biaxial stress field arises in them, which is summed 
up with the fields of residual welding stresses. Therefore, when selecting a technology for manufacturing critical 
welded structures, the results obtained during conventional uniaxial tests of samples are insufficient. The variety of 
factors affecting the performance of structures, and the difficulties of separate assessment of their influence, caused the 
need to maximize the approximation of experimental conditions to the real working conditions of the structure. 
Materials and Methods. Testing of full-scale structures has a number of advantages, but they are extremely expensive, 
and, as a rule, only one, the weakest link, is identified, the bearing capacity of the other structural elements remains 
unclear. For testing, UDI radiometric installations designed for different sample sizes were used. The presented 
installations allow testing samples of various shapes, types of welded joints (butt, T-bar), changing the position of 
welded parts. 

Results. Without rejecting the results obtained during the testing of full-scale structures in the works of the Bauman 
Moscow State Technical University, DSTU, NRC “Kurchatov Institute’ — CRISM “Prometey”, and the authors 
proposed to conduct the basic scope of the research on individual structural elements that would reflect the 
characteristic features of loading, manufacturing technology, and operating conditions. The design of the “fitting-sheet” 
connections was applied, which made it possible to increase the indicators of the failure initiation and propagation to the 
level of the base metal. 

Discussion and Conclusions. Schemes of structures for obtaining a biaxial tension or bending field in samples were 
presented. Samples tested according to the proposed schemes allowed us to draw conclusions about the performance of 
welded joints under conditions close to the actual operation of the structures under study. The proposed test scheme is 


used by research laboratories in our country and around the world. 


Keywords: tests of full-scale structures, design experience, installations for biaxial tension and bending, durability, 


“fitting-plate” connection. 
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Introduction. Thick-sheet hull pressure-operated structures are the most common type of high-duty welded 
structures. When these structures are loaded with internal pressure, a complex, predominantly biaxial stress state occurs, 
in which the fields of residual welding stresses are summed up with the stress fields from the external load. Under these 
conditions, there is often an increased sensitivity of both the base metal and welded joints to the presence of stress 
concentrators or a local change in the mechanical properties of the metal associated with the manufacturing process, 
which can cause a sharp decrease in strength. Currently, despite a number of fundamental studies on the strength of 
welded sheet and hull structures, there is still insufficient systematic data on the patterns of resistance to the initiation 
and development of destruction of large-sized elements of welded structures, depending on structural, technological, 
and operational factors. Therefore, for the choice of material and manufacturing technology of critical platework, the 
data obtained through conventional uniaxial tests of samples is sometimes insufficient. The variety of factors affecting 
the operability of structures and the difficulties of separate assessment of their influence necessitates bringing the course 
of the experiment to real conditions as close as possible. This has led to the spread of methods for testing full-size 
structures or their models. To conduct such tests in our country and abroad, stands equipped with sophisticated testing 
and diagnostic instruments have been created. One example is as follows: for testing full-size large-diameter pipes, 
there are stands in UralNITI, Volga Pipe Plant, VNIIST, CRISM “Prometey”, and in other laboratories. There are 
stands of a similar purpose in Japan. In the USA, there are more than 300 installations for testing structures operating 
under static and pulsating pressure. 

Tests of full-scale structures have a number of advantages: 

— real technology of manufacturing the structure is preserved, as a result of which the uncertainty associated with the 
influence of simulation of the manufacturing process on the test results is eliminated; 

—real scheme of loading of individual elements and components of structures is preserved; 

— they enable to establish the actual distribution of failures in the structure and to specify the accepted design 
scheme; 

— they provide assessing the resistance to the initiation and development of fracture of only one of the weakest 
structural elements; 

— they enable to study the effect of stress concentration and residual stresses on fatigue strength. 

Despite the undoubted value of field tests to assess the structural strength of the product, it is extremely irrational to 
use such tests to study individual factors affecting the strength of the structure. The abundance of factors simultaneously 
affecting the course of the study makes it difficult to analyze the causes of premature fracture and does not allow us to 
separately assess the degree of influence of each of them. As a result of testing of full-size structures, as a rule, it is 
possible to identify the weakest link, whereas the bearing capacity of the rest of the sections (elements) of the structure 
remains unknown. 

The above disadvantages can be practically eliminated if the bulk of the research is carried out on reduced-size 
models of shell structures. Such tests allow us to get answers to questions related to the structural strength and 
reliability of the product, to identify the degree of danger of the initiation and development of crack-like defects, to 
evaluate the effectiveness of various ways to increase durability. 

Without rejecting the test results obtained on simple uniaxial samples and when testing full-scale structures at the 
Bauman Moscow State Technical University, DSTU, NRC “Kurchatov Institute” — CRISM “Prometey”, the authors 
propose to transfer the bulk of the research to testing individual structural elements that reflect the characteristic 


features of its loading, manufacturing technology, and operating conditions. 
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In relation to shell and hull structures, during testing, it is required to take into account: 

— stress state (biaxial with equal or unequal components of the primary stresses); 

— loading condition (static or repeated static with different loading cycles); 

— environmental impact; 

— influence of operating temperature. 

Materials and Methods. This paper describes the experience of designing installations for testing metal and welded 
joints under biaxial stress. 

The analysis carried out in [1-5] has shown that the working conditions of metal and welded joints in steel-plate 
structures are most fully reproduced when tested under hydrostatic buckling. In this case, the sample is fixed or 
supported along the contour and loaded with hydrostatic pressure. 

The stress state that occurs in the sample metal depends on the shape of the sample, the conditions of its fixation 
(pinching or loose support on the die), the shape of the die. Figure 1 shows the loading schemes and the resulting biaxial 
bending or stretching stresses. 

When a flat sample supported along the contour of the circular hole of the die is loaded with hydrostatic pressure, a 
biaxial bend occurs. A significant part of the external convex surface of the sample is subject to uniform stretching with 
equal stress components 6; = 62 (Fig. | a). If the flat sample is securely pinched along the contour of the die hole, then 
biaxial stretching is applied to the biaxial bend. In the case of loading not a flat sample, but a spherical segment with a 
sufficiently large ratio of the diameter of the die hole to the thickness of the sample, the bending component is small, 
and it can be assumed that the central part of the sample is subject to the biaxial stretching with o2/o,=1 (Fig. 1 b). 
Biaxial stretching with unequal components within the ratio o,/o, =1.0...0.75 can be obtained by buckling the sample 


according to the scheme shown in Figure 1), using dies with elliptical holes. 
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Fig. 1. Loading schemes of samples to obtain a biaxial stress field in them : 
a — loading of a flat sample with equal stress components; b — loading of a spherical segment; 
c — loading of a sample in the form of a cylindrical panel 
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A further decrease in the ratio 6/0, = (0.7... 0.3) is obtained using the scheme shown in Figure | b, where a 
cylindrical-panel sample, pinched by a flange part between the cylindrical die and the cylindrical punch, is loaded with 
hydrostatic pressure [3, 4]. 

For steel-plate pressure-operated structures, two types of loading are characteristic: single (static) and low-cycle 
(repeated static). For the first, it is reasonable to use schemes a and b (Fig. 1). Under static loading, the scheme test can 
be carried out both on samples in the form of flat sheets and the pre-formed spherical segments [5—7]. 

Plates are preferable, since their production is less labor-intensive. Samples in the form of a spherical segment can 
reduce the edge effect of fixing the sample along the contour. However, the preparation of such samples requires plastic 
deformation, and this can cause a change in the mechanical properties of the material, which is not always correctable 
even by subsequent heat treatment. 

When testing under low-cycle loading conditions, all three schemes shown in Figure | can be used. However, 
preference should be given to biaxial bending according to scheme a, since it provides testing large thicknesses. 

The resistance to fracture nucleation and development under tension is usually evaluated using schemes 
b and c (Fig. 1) on pre-formed samples that are made in the form of a spherical segment or a cylindrical panel. 

The environment has a particularly strong impact on the results of long-term and repeated static tests. Of great 
interest is the resistance of materials under repeated static loading in corrosive environments [8—10]. 

To study the processes of nucleation and the development kinetics of fracture in the elements of hull structures, a 
series of installations has been developed in DSTU together with the NRC “Kurchatov Institute” — CRISM 
“Prometey”: UDI-550, UDI-980, UDI-1300. Their key characteristics for biaxial bending are presented in Table 1. 

Advantages of these installations: 

— powerful force actuators are not required for loading; 

—it is possible to test samples with a relatively uniform distribution of stresses on a large surface; 

— it is possible to test individual elements of welded structures with almost complete preservation of technological 


features and the pattern of stress distribution in the area of welds, to take into account the influence of the corrosive 


environment. 
Table | 
Characteristics of power units of installations for testing base metal and welded joints under biaxial bending conditions 
Geometry of test samples Yield steer put of Masmien pease. 
Installations di : ican materials MP 
iameter, mm ickness, mm 60.2, MPa a 
UDI-550 550 20-40 < 900 < 1,000 
UDI-980 980 30-60 < 900 < 1,000 
UDI-1300 1300 60-100 < 900 < 1,000 


The created installations make it possible to determine the cyclic strength of various materials in air and in a 


corrosive environment, to obtain information about the failure pattern with different structural and technological design 


of welded joints of full thickness, to evaluate the efficiency of various methods to increase durability. 


The design of the installations intended for testing under conditions of biaxial bending of welded joints made of 


sheet metal with a thickness of 20-100 mm with a yield strength up to 15 MPa is shown in Figure 2. They consist of 


four blocks 1, in which samples are fixed, a control panel 2 and pumping stations, which are not shown in the figure. 
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1,815 


1,400 


Fig. 2. Design of installations for testing under biaxial bending conditions: 
1 — test installation; 2 — control panel 


The basic design of the power units of UDI-550, UDI-980, and UDI-1300 installations is the same and is shown in 
Figure 3. 


—— 


Fig. 3. UDI-type power unit of installations for testing welded joints under biaxial bending: 


1 — die; 2 — hydraulic clamp; 3 — test sample; 4 — base; 5 — locking device (cup); 6 — inner ring groove; 
7 — ring locking sectors; 8 — protective casing; 9 — hydraulic jack 

The installation contains a die 1 and a hydraulic clamp 2, between which a sample 3 and a locking device 5 are 
placed, which is fixed on the base 4. The locking device is made in the form of a fixed smooth cup with an internal 
groove 6, in which annular locking sectors 7 (eight pcs.) are located along the entire perimeter, having an L-type shape 
in axial section. 

The installation works as follows. The sample 3 in the form of a round disk is mounted on the hydraulic clamp 2, a 
die 1 is mounted on top of it. Then, eight removable annular sectors 7 are laid in the groove of the locking device 5, 
pressing them tightly together so as to provide minimal gaps between them. 

The installation is covered with a protective casing 8. After that, the pumping station is turned on, and pressure is 
applied to the cavity A formed by the bottom of the locking device 5 and the hydraulic clamp 2. Under the impact of 
pressure, the hydraulic clamp 2 moves upwards, presses the sample 3 to the die 1, and the die — to the annular 
sectors 7. The sectors rest against the support surface of the annular groove and the inner surface of the cup 5, 
preventing the movement of the die 1, thereby fixing the sample 3 in the installation. 

Then, the liquid flow is directed into the cavity B under the sample 3. Under the influence of pressure, the sample, 
resting on the edge of the hole in the die 1, bends, and biaxial tensile stresses arise on its outer side. The magnitude of 
tensile stresses on the surface is regulated by the pressure in the cavity B. The hydraulic system allows for repeated 


static loading of samples with a loading frequency of ten cycles per minute. 
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The scheme of disassembly and extraction of the sample is shown in the right part of Figure 3. 

After the end of the tests, a hydraulic jack 9 is installed on the sample 3, which rests with its paws on the projections 
on the annular sectors 7. When connecting the jack to the hydraulic system of the installation, its piston presses on the 
sample and moves it together with the hydraulic clamp 2 down to the start position, which provides disassembling of 
the installation. 

In the process of testing, the number of cycles before the failure (crack length 8-12 mm), the crack development 
kinetics, and the destruction of the sample (in this case, it is loss of tightness) are recorded. 

The electrical and hydraulic circuits of the installations allow for static or repeated static loading. The loading cycle 
can be pulsating or with a positive coefficient of asymmetry. The frequency of loading samples is adjustable from two 
to ten cycles per minute. 

In cases where the effect of a corrosive medium is of interest, it is poured into the cavity formed by the sample and 
the die 1. To protect against corrosion, the die should be coated with a layer of epoxy resin or varnish, or an annular 
rubber collar should be glued to the sample, preventing the spreading of the corrosive medium beyond the test part of 
the sample. 

Depending on the testing purpose, samples of various types and sizes are used. Some of them are shown in Figure 4 


as an example. 


a) b) c) 


Fig. 4. Samples providing for assessment of the corrosive environment impact (3% NaCl solution): a — from the basic; 
b — linear notch; c — annular notch 


To assess the resistance to the nucleation and development of the base metal fracture in a biaxial stress field under 
the simultaneous action of a corrosive medium, it is reasonable to use smooth and notched samples. Comparative tests 
of these samples in air and in a 3 % NaCl solution have shown that this medium accelerates the processes of nucleation 
and development of fracture. When tested in air, as a rule, one crack nucleates and develops. In the presence of a 
corrosive medium (3 % NaCl solution), even in the presence of stress concentrators, a multifocal pattern of the 
nucleation and development of cracks is observed, which significantly accelerates the destruction process — 
by 2—2.5 times [11-13]. 

The paper considers a method for increasing low-cycle fatigue of welded joints of the “fitting — plate” type. It is 
proposed to increase durability by reducing stresses on the stretched surface, where a corrosive environment acts, 


through reducing the average stresses in the joint [14, 15]. Figure 5 shows the proposed welded joint. 
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Fig. 5. Section of the welded joint of “fitting — plate” type: 1 — fitting; 2 — plate (shell) 


The sample was tested on UDI-550 installation in a 3 % NaCl solution at rated voltages of 850 MPa. The test results 
showed that the number of cycles before the failure initiation doubled, and cracks appeared on the base metal. Photos of 
macrosections of the test results are shown in Figure 6. 


a) b) 
Fig. 6. Test results of “fitting — plate” welded joints in a corrosive environment: 
a — fitting connection with symmetrical welds; b — fitting connection 
with reinforcing weld on compressed fibers (photo of the authors) 


Analysis of the test samples allows us to establish that the durability of the fitting connection made by the proposed 
method turns out to be longer. In addition, due to the reinforcing influence of the weld metal from the side opposite to 
the action of the corrosive environment, the magnitude of the maximum tensile stresses in the most dangerous zone 
decreases; therefore, the resistance to the fracture initiation and development increases [16]. 

Research Results. Using the experience of designing individual installations and assemblies intended to test welded 
joints under biaxial bending conditions, scientists of DSTU and NRC “Kurchatov Institute” — CRISM “Prometey” have 
developed a series of installations designed to test sheet metal welded joints with a thickness of 20—100 mm with a yield 
strength up to 15 MPa under the biaxial bending conditions [17]. 
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Discussion and Conclusions. |. Equipment and testing methods of welded joints of the “fitting — plate” type, which 
operate under biaxial axisymmetric bending in a corrosive environment, were developed. This allowed us to evaluate 
the resistance to the fracture nucleation and development. The large size of the samples and the preservation of welding 
technology brought the results of these tests closer to the assessment of the bearing capacity of full-scale structures. 

2. To assess the feasibility of using new technological or design solutions, it is required to test samples of full-scale 
thickness. 

3. Failure of welded joints made according to serial (factory) technologies starts simultaneously at several sites 
along the fusion line, where there is a combination of unfavorable factors. The main role is played by factors of local 
stress increase, residual welding stresses, and the presence of corrosive environment. 

4. To increase durability, designs of welded joints with an asymmetrical arrangement of welds relative to the plane 
of the sheet are proposed, which have a smaller weld leg from the side of the action of tensile stresses and corrosive 
environment than from the opposite side. 

5. Installations and test methods have been implemented in the People's Republic of China and in the Republic of 
India. 
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Abstract 
Introduction. Determining the tumor (nodule) characteristics in terms of the shape, location, and type is an essential 


step after nodule detection in medical images for selecting the appropriate clinical intervention by radiologists. 
Computer-aided detection (CAD) systems efficiently succeeded in the nodule detection by 2D processing of computed 
tomography (CT)-scan lung images; however, the nodule (tumor) description in more detail is still a big challenge that 
faces these systems. 

Materials and Methods. In this paper, the 3D clustering is carried out on volumetric CT-scan images containing the 
nodule and its structures to describe the nodule progress through the consecutive slices of the lung in CT images. 
Results. This paper combines algorithms to cluster and define nodule’s features in 3D visualization. Applying some 3D 
functions to the objects, clustered using the K-means technique of CT lung images, provides a 3D visual exploration of 
the nodule shape and location. This study mainly focuses on clustering in 3D to discover complex information for a 
case missed in the radiologist’s report. In addition, the 3D-Density-based spatial clustering of applications with noise 
(DBSCAN) method and another 3D application (plotly) have been applied to evaluate the proposed system in this work. 
The proposed method has discovered a complicated case in data and automatically provides information about the 
nodule types (spherical, juxta-pleural, and pleural-tail). The algorithm is validated on the standard data consisting of the 
lung computed tomography scans with nodules greater and less than 3mm in size. 

Discussion and Conclusions. Based on the proposed model, it is possible to cluster lung nodules in volumetric CT scan 


and determine a set of characteristics such as the shape, location and type. 
Keywords: automated 3D Clustering, CT lung images, describing the nodule characteristics. 
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Introduction. Lung cancer is a leading cause of death worldwide among both men and women, with an impressive 
rate for nearly 10 million deaths in 2021. Around one-third of deaths from cancer are to smoking and second-hand 
smoke. Many cancers can be cured if detected early and treated effectively (1). 

Computed Tomography (CT) is one of the most commonly medical imaging techniques used in the screening of 
lungs by generating three-dimensional imaging modalities and shows the lesions that cannot be visualized by 
conventional chest X-ray. CT technology is used widely by clinicians to detect, analyze, and diagnose numerous 


asymptomatic disease lung diseases such as pulmonary nodules and lung cancers, that cannot be detected by other 
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medical imaging technologies. Where, early detection of these diseases is an important of health promotion and chronic 
disease prevention. However, it is noted that a big challenge associated to radiologists is the analysis of huge amount of 
data generated by CT technology. Therefore, computer-aided diagnostic (CAD) is needed to provide a computerized 
diagnosis to assist clinicians for detection and interpretation of CT lung scans (2, 3). Segmentation and clustering are 
tasks that manipulate the medical image to extract significate information and reduce the search area of the region of 
interest in image components. Where, using clustering on volumetric scan of CT can provide more details about nodule 
shapes, characteristics, and their types which are suspected be missed by radiologists during the diagnosis. Additionally, 
it provides a precise description of lung’s tissues with low-intensity (4, 5). The outcomes of 3D clustering can improve 
the precision in diagnosis of lung by identifying the nodule structures in the different levels, and the number of 
nodules (6, 7). This study aims to improve the accuracy of three- dimensional clustering of lung’s nodules by using 
volumetric CT scans of lung. Where, the three-dimensional clustering efficiently can be used to observe the current 
nodule’s progress which serves significantly the early diagnosis of lung. Furthermore, the precise details that can be 
extracted from the three-dimensional clustering offers a high potential to detect nodules and lesions with less 
than 3 mm. 

Many studies focused on the possibility of developing a computer-aided system for diagnosing, detecting, and 
segmenting lung lesions in CT scans of lung (8). However, without using volumetric CT scans of lungs, the accurate 
identification and characterization of small pulmonary nodules have not been defined in the studies. Where, 
implementing quantitative analysis of CT scans of lungs face an obstacle due to inter-scan image intensity variations 
and irregular shapes of lung lesions. A major difficulty is encountered during clustering of lung’s lesions because of the 
complex characteristics of lung’s lesions and overlapping in intensity distributions in such cases. Additionally, some 
parts of lung’s lesions may not be distinguished correctly because of the image resolution and the complexity of 
diagnosing lung’s lesions in a visual examination. Therefore, a multi-dimensional CT scan is required to determine type 
and volume of lesions. 

Javaid, Javid (9) proposed an automated method for segmenting and recognizing the pulmonary nodules into juxta- 
vascular or juxta-pleural. Where, thresholding and k-means were combined to segment lung and detect the pulmonary 
nodules respectively. Then, the pulmonary nodules were classified into sex groups according to percentage of 
connectivity with lung walls and thickness. Fetita, Preteux (10) developed a method for discriminating lung’s nodules 
from other dense structures by using volumetric space of the thorax. Where, the achieved false positive rate was 8.5 per 
exam. Fetita, Preteux (10) applied a three-dimensional active contour method to segment 96 lung’s nodules, then the 
texture features were extracted and classified into classes by a linear discriminant analysis (LDA). Ozekes and 
Osman (11) suggested a CAD system based on volumetric CT scans of lung to detect lung’s nodules which have 
diameters (3.5—7.3) mm. In this study, four classifiers were used to evaluated the suggested system. 

Furthermore, El-Baz, Elnakib (12) exploited the efficacy of genetic algorithm (GA) to isolate nodules, arteries, 
veins, bronchi and bronchioles from the surrounding anatomical structures based on volumetric and two-dimensional 
CT lung scans. It was proved that the lung nodules were recognized precisely when using volumetric CT scans. In this 
study, it was noted that the false positive rate was reduced significantly when utilizing volumetric CT scan of lung 
compared to two-dimensional CT scan. The achieved sensitivity and false positive rates were 82.3 % and 9.2 % 
respectively. 

In this study, we proposed a new system to describe automatically lung’s lesions in volumetric CT scans, as well as, 
determining the nodule numbers, type location, shape, progress of lesion. The main contribution of this study is to 
discriminant the complicated lung’s nodules that have non seen through the clinical routine. A successful system can 
improve the accuracy of the diagnosis process. The following research objectives focus on how the above aim will be 
achieved: 

— to extract the most distinct characteristics of the lung’s nodule to that helps to determine treatment decisions and 
predict expected outcomes; 

— to detect the complicated lesions cases that cannot be identified with certainty by radiologists; 


— constructing an efficient CAD system that can determine precisely lung’s nodules number and types. 
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The rest of this paper is organized as follows. In Section 2, material and method will be described and the 
conclusion will be discussed in Section 3. 

Materials and Methods. The core-aim of the present study was to build an automated algorithm of clustering 
nodules from CT lung scans in an axial plane. The proposed algorithm is divided into three stages as shown in Fig. 1. 
First, the CT lung scan undergoes an enhancing process and separating of lung from background in the pre-processing 
stage. Then, the CT lung scan is clustered in 2D and 3D for detecting the vessel and nodules from normal tissue in the 
post-processing stage, as well as reducing the area region about nodules in 2D clustering while describes the nodule 
shape in 3D clustering. The results of 3D clustering are evaluated using two ways (3D-DBSCAN method and Plot.ly 
application). The evaluation stage confirms the outcomes precision of previous stage (post-processing) in detecting the 


number and shape of nodules that founded in the complex case, which is recorded in radiologist’s report. 
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Fig. 1. The Proposed Model of Lung CT scan Clustering 


Data Collection This study included 85 lung CT scans with pulmonary nodules, were downloaded from the TCIA 
(Cancer Imaging Archive) website (13). The provided data was annotated and screened by four radiologists to 
determine and extract the diagnostic details of lung’s nodules such as number, size, location, volume, diameter, type of 
lung’s nodules. In addition to slice number of the pathological lung. Figure 2, shows samples of CT lung scans form the 


provided dataset in axial viewing. 
_w=t 


J 


a) b) 
Fig. 2. Lung CT slices in the axial view: a — healthy CT lung scan; b — pathological CT lung scan 


Pre-processing of Lung CT Scan. To reconstruct the CT images, a million of independent detectors are used to 
collect the x-ray signals. Subsequently, the CT artifacts may be occurred during the reconstruction process. Where, 
these artifacts cause intensity variations across consecutive re-constructed CT scan slices. Therefore, the lung CT scans 
are pre-processed by implementing a set of algorithms such as image enhancement by gaussian filter, intensity 
normalization by histogram normalization method, lung extraction by thresholding (1, 14) as shown in Fig. 3. 
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Fig. 3. Lung boundary identification samples: a — Original CT image; b — Preprocessed CT image 


Three-dimensional Clustering. In this step includes implementation of a set of image pre-processing methods to 
prepare the CT scans of lung more appropriate for implementing the K-mean clustering method. Where, reduction of 
parenchyma tissue was achieved by isolating the intensity values that represents the vessels and nodules together. Then, 
vessels and lung’s nodules were separated by implementing the K-mean clustering method on each CT image 
separately. It is an efficient and non-supervised technique that allows the cluster center to shift and fit the ROI (15-17). 
Subsequently, any hole appearing in the clustered CT scan of lung is removed by implementing a set of morphological 
operations (1). Then, the clustered CT images of each patient were stacked into single a volumetric CT scan and 
underwent the 3D graphics functions later. 

Finally, a three-dimensional volume reconstruction from cross-sectional CT scans is determined by using isosurface 
function which was provided by MATLAB and used to measure the geometric computations of clustered data, and 
improved volumetric visualization of objects. Where, pixel coordination (x,y,z) and intensity values of pixel of each 
cross-sectional clustered scan was given as input to the function. Moreover, isocaps and isonormal functions were used 
to enhance produced the three-dimensional representation of the lung’s nodules and vessels. Consequently, the three- 
dimensional clustering stage detects nodules, shape, type, and progress and finds out the unseen lesions in the CT scan 
of lung. The suggested algorithm was applied to cluster nodules and vessels in volumetric CT scan of lung. 

First step was the threshold, where each value of the intensity in the volumetric CT scans by a specific value to 


extract correctly the lung’s nodules from the CT scan of lungs by using Eq. 1. 
1 ifP@,y) >T 
Gy) = { 0 Peden (1) 

Where, P(x,y) indicates the intensity values representing the nodules and vessels which are greater than the 
threshold value which is set experimentally. 

The second step used K-means clustering on each CT lung images individually to eliminate any undesirable holes 
that occurred in the clustering process. The output of clustering process was labelling pixels of CT images into two 
main clusters (nodules and vessels). K-means clustering implemented by using the following steps. 

1. Initialize k to specify cluster center, C=c,, C2,..Cx. 

2. Determining the distance between CT lung image’s pixels and the cluster centers by using Eq. 2: 

S = Dis Dyas Uy — gy UP (2) 
where, x; is an intensity pixel associated with m and n coordinates. 

3. Assigning the CT lung’s pixel with minimum distance to the nearest cluster’s center. 

4. Updating the cluster centers according to Equation (3): 

Ci = a Xij (3) 
where, c; and N;, ; indicate the i“cluster center, and all pixels of the i” cluster respectively. While, x;,,; signifies all 


pixels of i” cluster center. 

5. Repeating steps 2 to 4 until minimizing the objective function S to the lowest value. 
Finally, implemented a three-dimensional function on the labelled clusters of CT images to generate the volumetric 
data. Where, the volumetric data was plotted according to a set of calculations. The following procedure was used to 
combine the labeled clusters and produce lung’s nodule’s shape and type (18). Algorithm 1 shows the pseudo-code for 


the generating of volumetric lung’s nodule steps. 
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Algorithm 1: Pseudo-code for generating of volumetric lung’s nodule. 


Input: Clustered image (x.y.z). 
Output: Output volumetric figure (i,j,k). 
Begin 
Apply Jsosurface function to present a volume visualization. 
Apply Jsocaps function to adjust colors of faces of the produced volumetric lung’s nodules. 
Apply Azimuth and Elevation functions to create axes in three-dimensional view. 
Apply Zsonormals to determine the normal vertices and create a surface shadow through 
vertices gradient. 


Results. The plotted three-dimensional space of lung’s nodules enhanced visualization of nodules and vessels and 
provided more descriptions on the lung’s nodules than what extracted by two-dimensional space. Where, the former 
provides more details regarding shape and locations of lung’s nodules. Figure 4, shows an example of pathological lung 
that has two nodules with a red referral that are appear between slices no. 60 and 64. After implementing the three- 
dimensional clustering, a thin tissue between lung’s nodules that is invisible by two-dimensional clustering, becomes 
visible to the radiologists. Furthermore, it extends through CT slices no. 61 to 63. Therefore, by depending on the two- 
dimensional clustering leads to make the inspection process by the radiologist prone to error and may diagnose one 
lung’s nodule as two separated nodules. 

Moreover, the proposed system provides more details about the lung’s nodule type. Where, each lung’s nodule was 
classified into juxta-vascular and juxta-pleural nodules. Figure 5, reveals how an early lung’s nodule was detected and 
notified and this will help the clinicians to improve patient survival rates. Figure 6 shows how the lung’s nodules are 
connected through their center pixels. 


ee OT EE 


e) dD 


Fig. 4. Original CT Lung scan between Slice 60 and 64 with a Red arrow to indicate Two Attached Nodules: a — Candidate Nodule 
in the 60th slice; b, c, d— lung Nodules through slices 61, 62, 63; e — Nodule in slice 64; f— Normal slice after slice 65 
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Nodule in 


~ 60 slice 


Fig. 5. Three-dimensional Clustering of CT lungs Attached Nodules in 60th and 64th Slices 


350 #06 


Fig. 6. Three-dimensional Clustering of two lung nodules that are located in different slices 


Although the three-dimensional graphic is evaluated visually, it is essential to evaluate quantitively the quality of the 
achieved results. The Density-based spatial clustering of applications with noise (DBSCAN) method and 3D plot.ly 
application were used to evaluate the performance of the proposed three-dimensional clustering method (19, 20). 

The DBSCAN method is one of the most popular and standard clustering techniques (21—23). It was used to set of 
nearby pixels are grouped to gather into a single cluster, and isolate the outlier pixels in a low-density area. The 
clustering process of this method depends on two key parameters (Epsilon and MinPts) that deal with core point in the 
volumetric space. Where, Epsilon is the distance that specifies the neighborhoods. Two points are considered to be 
neighbors if the distance between them is less than or equal to eps. While, MinPts represents minimum number of data 
points to define a cluster. Then, a point is considered as a core point if there are at least MinPts number of points 
(including the point itself) in its surrounding area with radius eps. 

DBSCAN method was modified to cluster the volumetric CT lung scans to compare its output with the proposed 
system in this study by comparing the cluster’s centers that were achieved by both methods. Where, DBSCAN inspects 


and connects all objects to the core pixels based on eps and MinPts measures. Synchronously, all the points that do not 
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represent core points are ignored. Many advantages are achieved by using this method, first, no need to initialize cluster 
numbers like k-means cluster method. Secondly, random shaped clusters are discovered easily. Finally, if the data is 
fully understood, it is easily to specify the MinPts and ets parameters. 

Fig. 7 shows how the lung’s nodules were linked through CT slices. Where, the 3D-DBSCAN method evaluated the 
connection between attached nodules compare with the clustered lung’s nodules of the proposed system in this study. It 
is noted that the center’s pixels were approximately matched as shown in Fig. 8. Fig. 9, shows how the lung’s nodules 


were linked and spread by a thin line between CT slices no. 60" and 64". 
Noise 
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Fig. 7. 3D-DBSCAN clustering Method, yellow cluster refers to Linked Nodules, 
and other colors refer to the Remaining Nodules of the Same Patient 
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Fig. 9. DBSCAN Method Just Clusters the Linked Nodules and Proves One Nodule 


Another 3D application was utilized in this study to evaluate the proposed system. It is named plot.ly application. 
This application used colored bubbles to locate pixels within the candidate nodule and its structures in different slices, 
as shown in Fig. 10 and Fig. 11. It was initialized by the coordinates of lung’s nodules centers and it extended 
consecutively through similar structures in neighbor CT slices. Moreover, it was essentially to initialize pixel size and 
pixel spacing of CT scan, which are determined by Eq. 4 and Eq. 5. 


: i =o = _ ((Field of view(FOV)) 7 _ (360mm 2 = 2 
Pixel spacing(x — size, y — size) ( (aeesne ) ( a ) 0.703 (4) 
voxel size = pixel spacing (x — size) x pixel spacing(y — size) x thickness = 0.703 x 0.703 X 2.5 = 1.235 (5) 


Where x and y indicate the distance from the center of one pixel to the center of an adjacent pixel in the X and Y 
axes respectively. It is determined by squared ratio of the CT image field of view to the image array size (512X512). 
The image parameters of the provided dataset; pixel spacing, thickness, and voxel size were (0.703/0.703), 2.5 and 
1.235 respectively. Fig. 12, shows the centered pixels of lung’s nodules were labeled with two colors: red and green that 
represented actual center that was recognized by radiologist and the determined center by the proposed system. The 
black and blue colors show the overlapped area between clusters in different slices, and the remaining pixels of the 
cluster areas respectively. Generally, the plotted volumetric lung’s nodules by the Polt.ly application and 3D-DBSCAN 
methods proved the efficacy of the suggested system in this study. Where, all used methods described accurately the 
connected regions between linked nodules in CT lung scans. 
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Fig. 10. Three-dimensional Plot of lung’s Nodules was distributed on Different Levels from Slice 60 to 64 
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Fig. 11. Three-dimensional Plot of lung’s nodule of the same case in Fig.10 but from different view, Black Bubbles Represent the 
Overlapping between Nodules to Show the Connected Regions 
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Fig. 12. The Overlapping between Centre Pixels of Clustered and Actual Nodule in the 60 Slice with 
Two Colors Red and Green Respectively 

Discussion and Conclusions. The main objective of this paper was to build and implement a fully automated 
algorithm to describe the pulmonary nodules in CT lung images regarding shape and location matching the radiologist’s 
report. The proposed method could cluster the nodules and showed their shape, location, type, and progress within the 
lung organ. The final algorithm outcomes and case study were evaluated by other 3D clustering methods emphasizing 
the accurate research results. This work is performed within several stages, with different methods are used for this 
purpose. Firstly, threshold and K-means clustering segment the nodules and vessels of CT lung images as objects with 
high intensity in 2D images. With 3D function, the 2D clustered images are stacked to display the nodule progress in 
3D view. The algorithm results were validated by matching the center pixels of detected nodules by 3D Clustering and 
actual nodules reported by radiologists. In addition, the study identified a complicated case in the data, which is 
suspected the radiologists may miss because it was invisible. In this case, the algorithm proves and shows an attachment 
between two nodules making the nodules as one nodule against the radiologist’s report. The evaluation stage has been 
done by applying the 3D-DBSCAN and 3D application methods, which confirm the proposed method results precision 
for this work. The Clustering of 3D-DBSCAN demonstrates the number and location of nodules and reflects the nodule 
shape as far as to describe the type of nodule. Also, determining this case confirms the work abilities in the early 
detection of nodules, which provides a chance for recovery patient. Our algorithm has a promising result on the 
standard data of CT scan images in detecting the number and location of nodules, which correspond to radiologist 
report, and automatically providing nodule shape and type. 
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Abstract 
Introduction. One of the most significant tasks facing road experts is to maintain the transport network in good 


condition. The process of selecting an appropriate approach to providing such condition is quite complex since it 
requires considering many parameters, such as the existing condition of the pavement, road category, weather 
conditions, traffic volume, etc. Recently, the rising trend of digitization in the industry has contributed to the use of 
artificial intelligence to address problems in several fields, including the bodies in charge of operational control over the 
status of roadways. Within the context of any control system, the main task of the control system is to carry out reliable 
forecasting of the operational state of the road in the medium and long term. 

Materials and Methods. This study investigated the possibility of using artificial neural networks to assess existing 
pavement characteristics and their potential application in developing road maintenance strategies. A back-propagation 
neural network was implemented, trained using data from 1,614 investigated sections of the M4 “DON” highway in the 
road network of the Russian Federation in the period from 2014 to 2018. Several models were developed and trained 
using the MATLAB application, each with a different number of neurons in the hidden layers. 

Results. The results of the models showed a convergence between the inferred paving state values and the actual values, 
as the multiple correlation coefficient (R’) values exceeded 92 % for most of the models during all learning stages. 
Discussion and Conclusions. The findings suggest that public road authorities may utilize the established models to 
choose the best road maintenance strategy and assign the most efficient steps to restore road bearing capacity and 


operation. 
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maintenance, pavement management system. 
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Introduction. Road infrastructure is one of the resources that are vital to any community, so it is necessary to think 
about ways to maintain its efficiency. However, roads deteriorate over time as a result of being subjected to various 
deterioration mechanisms [1]. As a result, it is essential to manage these resources strategically to maximize their 
service life and technical indicators while also lowering maintenance costs [2]. 

Pavement management systems (PMS) aid in the maintenance of pavement performance at an appropriate level for 
use and provide a variety of financial and social advantages. Moreover, it contributes significantly to the cleanliness of 
the surroundings during the usage or maintenance [3]. Accordingly, if the condition of the pavement network is allowed 
to deteriorate, the damages will not only increase the costs of road maintenance to raise their efficiency. But, it will 
cause environmental harm due to emissions and traffic noise, as well as expose travelers to the hazards of automobiles 
travelling on rough roads [4]. 

Due to the significant damage caused by the deterioration of roads, transportation agencies in various countries 
allocate a large part of their annual budget to the maintenance or rehabilitation of these damaged roads. The importance 
of maintenance planning work is to choose the appropriate time and method to reduce the rate of deterioration of 
pavement sections according to the available budget. To reach this goal, it is required to develop models that help 
improve the development process and take into account the factors affecting the maintenance decision making [5]. 

Multi-criteria models are commonly used in complicated decision-making situations, like those involving the 
network-wide maintenance of road infrastructure [6]. Multi-criteria models emerge as a scientific approach for 
evaluating multiple alternative interventions in the framework of road rehabilitation, taking into account the properties 
of pavement sections [7]. 

It is well recognized that a pavement's surfacing behavior reflects its condition; for this reason, pavement 
administrations depend on knowing the value of the PCI index through the process of visual inspection of the roads 
periodically in order to extend the life of the pavement and improve service levels [8]. The value of the pavement 
condition coefficient is one of the most significant factors based on which the most appropriate method of maintenance 
is determined, in addition to other factors, such as weather factors and traffic volume. 

Over the last two decades, it has become clear that traditional knowledge of road maintenance is not enough, and 
there is an urgent need to develop new methods that help collect and process data for use in pavement management [9]. 
Due to the rapid development in the domains of information technology and artificial intelligence networks, many 
opportunities have been provided to find solutions to many problems in various fields, including road engineering. 

Artificial neural network models are an application of the field of artificial intelligence used to solve non-linear 
geometric models, such as forecasting, recognition, and estimation of different patterns [10]. Artificial neural network 
models mimic the human brain's ability to solve issues using previous experiences. As a result, we can use this 
approach for choosing the most appropriate way to maintain the road after calculating its current or future condition 
instead of relying on human experience and knowledge. 

Artificial neural networks (ANN) review: 

Many studies focused on the possibility of using the artificial neural networks approach to determine the current and 
future state of the pavement to using it in the development of pavement management systems and determining the 
method of maintenance. We will mention some of these studies as follows: 

R. Kumar, et al. created a neural network model for determining the condition rating of flexible pavements 
depending on pavement distress data collected. The study result appeared that the neural network was able to calculate 
pavement condition rating accurately [11]. 

Hamdi, et al. built an artificial neural network model using distress data generated from a visual assessment method 
to determine the Surface Distress Index. The positive modeling results showed that it could be used to determine the 
SDI coefficient [12]. 

After the results of the ANN models showed superiority in predicting the values of the pavement condition, many 
researchers were interested in developing models that can predict the values of pavement distress as an indicator of the 
pavement condition. 

D. T. Thube used a neural networks approach to construct four models to predict pavement defects (cracking, rut 
depth, roughness, and raveling) as pavement condition indicators. The small differences between the observed and 
calculated distress values were evidence of the model's success [13]. 

L. Yao, et al. used a neural network approach to construct five models that predict five indicators of pavement 
distress: the transverse crack index, the skid-resistance index, roughness, the rut index, and pavement surface distress. 
The results of the models showed an acceptable performance in predicting the condition of the pavement, with an 
average testing R” of 0.8692 [14]. 

A. Shtayat, et al. [15], M. Mazari and D. D. Rodriguez [16] studied the possibility of predicting the International 
Roughness Index (IRI) based on distress data or different pavement performance indicators, such as traffic coefficients 
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and structural properties of pavement through neural network models or combined with gene expression programs or 
regression analysis. However, the results of neural network models showed an acceptable superiority compared to other 
methods. 

A.M. Mosa suggested a neural network-based approach for diagnosing pavement distresses and optimizing 
solutions for maintenance techniques suggestions. The created system supplied the optimal solutions taking into 
consideration technological, economic, and environmental requirements [17]. 

The results of the mentioned literature review showed the ability of neural network models to represent the complex 
and nonlinear relationship between several variables. Therefore, it can contribute to solving road engineering problems. 
We also note the use of most previous studies of ANN models with an inverse propagation logarithm, and we will 
mention a summary of them in the following section. 

Description of back-propagation artificial neural networks: 

The back-propagation algorithm is one of the most widely utilized algorithms of artificial neural network models. It 
has a high capacity and speed in processing new data after conducting several iterations for a set of data during the 
education stage [18, 19]. 

The structure of the back-propagation algorithm for artificial neural networks consists of three or more layers, 
including the input layer, the output layer, and one or more hidden layers, as shown in Figure 1. The function of the 
neurons inside the input layer is to receive data and then pass it through the hidden layers, which have the ability to 
examine the correlations between variables and process the data to the output layer, which provides the results by the 


neurons in it [20] . 
Inputs 


Input layer 


Neuron Hidden layer 


Connection 
Output layer 


Outputs 


Fig. 1. Artificial Neural Network Construction Layers 


In order for the neural network to perform the function assigned to it, the data must pass through three stages. The 
first is the process of training the network, in which it understands the effect of each variable and its relationship to 
other variables. Then we move on to the second stage, which is the validating process of the model, and finally, the 
network testing stage [21]. 

The back-propagation approach employs supervised learning, which implies that we supply examples data for the 
inputs and outputs to the network for learning [22]. The model assumes random weights for each variable and then 
passes the data forward to calculate the output. After calculating the difference between the computed output and the 
target values during the back pass, the weights are adjusted again to reduce the error rate and return the data from the 
output layer through the hidden layer to the input layer. The described learning process is repeated several times until 
the learning process is stopped when the error reaches an allowable value [19]. 

The weights determined during the training phase are employed in the testing process to calculate the output of a 
new input data set that was not used during the learning process. An evaluation of the ANN model is based on the 
amount of computed error. 

The research scopes: 

The study aims to use the field of artificial intelligence to build a model that helps engineers or decision-makers who 
prepare a road maintenance map to know the condition of the pavement without the need to conduct a periodic visual 
inspection. This reduces the dangers to individuals in charge of the examination process, as well as the effort, time, 
expense associated with it, and the errors that occur throughout the measuring procedure. To achieve this goal, we 
utilized the values of the deflection measurements coming from the falling weigh deflection test and the pavement 
condition index values from the Russian Road Company project for developing the M4 highway between the periods 
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from 2014 to 2018. We prepared a database for training and testing artificial neural network models by MATLAB 
software. The importance of this research is to accurately describe the future condition of the asphalt pavement to use 
the appropriate maintenance method to slow down the rate of asphalt deterioration. 

Materials and Methods 

Preparation of the database: 

The data obtained for constructing the ANN model were collected from a study carried out by the Russian Road 
Corporation for developing the M4 highway between Moscow and Krasnodar. The available data includes inspection 
date, air temperature, asphalt layer thickness, base layer thickness, traffic volume, and precipitation rate and deflection 
values based on the falling weight deflectometer test. 

The database was used to calculate three criteria that express the condition of the asphalt surface, which are the 
probability of success of the pavement sectors against the following three defects: fatigue, roughness, and rutting, as 
well as the modulus of elasticity of the pavement layers (asphalt layer, base layer, and subgrade layer) that were 
determined using a software package PRIMAX as input data and Regarding calculating the input parameter of 
“pavement surface life”, it was determined by taking the difference between the deflection survey date and the date of 
construction or last rehabilitation date. Also, the initial pavement condition mentioned in the same study, represented by 
PCI values, was used as output data to train ANN models to predict PCI values. 

Calculating pavement defects using the logistic model: 

Non-destructive testing is an important approach for evaluating pavement structures and is widely accepted as a 
reliable way of determining the structural condition of existing pavements [23]. The falling weight deflectometer 
(FWD) is well-known for its effectiveness in determining the structural condition of pavement and assisting in defining 
the best treatment option possible, which reduces the deterioration of the pavement [24, 25]. 

The values expressing the resistance of the pavement sections to the above three pavement defects were calculated 
using the logistic model equations developed by the Federal Highway Administration Long-Term Performance (LTPP). 

Equations 1 and 2 show the general formula, the linear formula, that were produced from the logistic model to 
calculate the three distress indicators of the pavement for all sectors used in building the ANN model. 


Equation (1) described the general formula of the logistic model. 
/ 


p(event)= —— (1) 
Equation (2) described the calculation formula of the exponent term in the general linear. 
b agt a,b; t ayb Pel 4,0, (2) 


where p represents the probability of an event occurring, b represents the exponent variable and (a, a), a2. a,) are the 
constants of the variables of the linear equation. 


The variables that influence pavement fatigue to determine the value of the parameter (b) are: 

— D, denotes the deflection measured in the loading plate's midpoint; 

— AADTT; 

— the base layer type. 

To compute the exponential component in equation (1), we utilized the variable values generated from the logistic 
model for pavement fatigue cracking in Table 1. 


Table 1 
The coefficients values utilized in the linear equation to calculate the condition of fatigue cracking 
Variables (bi) aj 
T 154.764 a, = coefficient values of the 
1 : . 
AADIT —0.0005073 eee 
I, = The utilized deflection value's 
Pavement type 0.3774 : 
mutual index 
Constant —0.2202 
Equation (3) is used to find the variable J; [26]. 
1 
lim> 3) 


The variables that influence pavement roughness to determine the value of the parameter (b) are: 

— D, is the measured deflection value at distance 200 mm from the center of the loading plate; 

— the volume of trucks in the Class 9 classification; 

— the pavement surface age from construction or from the latest maintenance. 

To compute the exponential component in equation (1), we utilized the variable values generated from the logistic 
model for pavement roughness in Table 2. 
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Table 2 
The coefficients values utilized in the linear equation to calculate the condition of roughness distress 
Variables (bi) aj 
I 739.849 a; = coefficient values of the equation. 
2 : a : j 
Guncnche 9.189 lh= _ utilized deflection value's 
mutual index. 
Class 9 volume —0.0006781 
Constant 0.8375 
Equation (4) is used to find the variable J, [26]. 
1 
b= — (4) 


The variables that influence pavement rutting to determine the value of the parameter (b) are: 

— D3 is the recorded deflection at a distance of 300 mm from the load plate's center; 

— D4 is the deflection recorded at 450 mm from the load plate's center; 

— the volume of trucks in the Class 9 classification; 

— average annual precipitation in the zone area (mm). 

To compute the exponential component in equation (1), we utilized the variable values generated from the logistic 


model for pavement rutting in Table 3. 


Table 3 

The coefficients values utilized in the linear equation to calculate the condition of rutting distress 

Variables (bi) a; 
Cl; —0.01146 aj=coefficient values of the 
Precipitation (mm) —0.0005259 equation. 
Class 9 volume —0.0007688 CI; = Curvature index for the i” 
Constant 2.6586 deviations values. 
Equation (5) is used to find the variable CT; [26]. 

Cl;= D; = Dd, (5) 


Initial Database for Neural Networks 

The database for entry in the neural network is comprised of a set of input data and a set of output data. The input 
data are the measured values of the following parameters: pavement surface life, asphalt layer thickness, base layer 
thickness, p(fatigue), p(roughness), p(rutting), the elastic modulus of (surface layer, base layer, and subgrade layer) 
shown as their mean value for an individual pavement segment. The model outputs are the pavement condition index 
values for the pavement sections that are available in the database. The goal of this database is the preparation and 
sorting of data in a format appropriate for entering it into the neural network. 

Application of ANN 

For the objective of this study, artificial neural networks with back-propagation algorithms and different numbers of 
neurons (8, 9, 10, 11, 12, 13, and 14) in two hidden layers were applied using the MATLAB software. The total 
quantity of inspection data utilized in the construction of the ANN model was 1,614 representing 51 different sectors 
along the M4 highway. The data used in the training process included the minimum and maximum values of the 
variables to increase the efficiency of the models. 

The database was divided into three groups: a large amount of database in which the neural network conducts the 
learning process on 70 % of the randomly selected data; the second group represents the validation process with 15 % 
of the database; the remaining 15 % of the data for the third group was for testing the developed network. Figure 3 


shows the configuration for one of the created neural networks during training. 
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Fig. 2. The Applied ANN Model Structure 


Each of the input parameters is represented by one of the nine neurons in the input layer. Eleven neurons are applied 
in the hidden layers, split into seven neurons in the first hidden layer and four neurons in the second layer. The output 
layer includes one neuron that represents the number of output data. The previous figure shows the direction of data 
transmission from the input layer to the output layer through two hidden layers and vice versa to correct the error rate in 
the output values until the permissible limit is reached. 

Several numbers of layers and neurons inside the hidden layers were applied to increase the performance of the 
developed network. The training procedure was repeated several times until the optimal model that accurately expressed 
the relationship between the input and output data was found. The effectiveness of the trained models were compared 
using statistical analysis of the outputs, which were represented by the mean absolute error (MAB), and coefficient of 
multiple correlations (R’), Root Mean Square Error (RMSE), and mean absolute percentage deviation (MARD). 

Results. After the learning step, the neural network was subjected to testing. In the testing process, weights were 
fixed to values adopted at the end of the learning process. A new group of input and output data was offered to the 
network. Network output results were compared with required output and statistically analyzed. 

The following regression and performance plots illustrate the output of the network models concerning training, 
validation, and test sets. Throw the training process, the optimal values of the network elements were carried out after 
1,000 repetitions with (MSE = 0.0224). The best validation performance of the demonstrated model was 0.0344, which 
was achieved after 22 repetitions for the model with structure (9—10-1), as shown in Figure 3 (a, b). The optimal values 
of the network elements were carried out after 1,000 repetitions with (MSE = 0.0234). The best validation performance 
of the demonstrated model was 0.03171, which was achieved after 11 repetitions for the model with structure (9—11-1), 
as shown in Figure 4 (a, b). The low value of MSE of the selected models indicates that pavement conditions of roads in 
the same domain may be anticipated more accurately. 
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Fig. 3. Results of training an ANN model with a structure (9—10—1): a — regression graphs; b — performance chart 
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Fig. 4. Results of training an ANN model with a structure (9—11—1): a — regression graphs; b — performance chart 


According to the output results, notice a high correlation coefficient between the outcomes of the neural network 
computation and the needed outputs. There is a significant correlation between them where most of the data falls near a 
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45-degree line. As shown in the previous regression charts displayed for all stages, whether training, validation or 
testing phases, the correlation coefficients (R) values are greater than 95% for all stages, indicating the quality of the 
generated models and the ability to estimate the outputs correctly. 

Evaluation of Artificial Neural Network: 

The created ANN model may be saved as a MATLAB file once the training procedure is completed, and the desired 
accuracy is achieved. We can use it to do a forward computation and forecast pavement conditions for maintenance and 
rehabilitation purposes. 

For this purpose, a new dataset of 37 sections was assigned to asphalt pavement to evaluate the ability of the 
developed neural networks used in the previous section to predict the output values. In this step the models were 
defined by input data for the specified parameter values, which were: (asphalt layer thickness, base layer thickness, 
surface life, p (fatigue), p (roughness), p (rutting), E (asphalt), E (base), and E (subgrade)). As a result, the created 
networks were able to anticipate output data (pavement condition index- PCI) based on early experience. 

The criteria for choosing the best iteration for every model were: mean absolute error (MAE), coefficient of 
multiple correlations (R’), Root Mean Square Error (RMSE), and mean absolute percentage deviation (MARD) of 
the estimates (PCI). They ensured that the predicted values were within reasonable data limits, as shown in Table 4. 
Figure 6 (a, b, c, d, e, and f) presents the relationship between the prediction results of the best iteration for every model 
with different numbers of neurons in the hidden layers, and the target output results. 


Table 4 
Statistical assessment of the created models’ output (forecasted results) 
Statistical comparison Pavement Condition Index — PCI 
criteria 8n on 10n lin 12n 13n 14n 
MAE 0.1415 0.1242 0.1185 0.1317 0.1362 0.1594 0.1190 
R° 0.9204 | 0.9493 0.9451 0.9365 0.9342 0.8686 0.9483 
RMSE 0.2058 0.1651 0.1710 0.1840 0.1862 0.2765 0.1666 
MARE 3.4478 3.0430 2.8738 3.1545 3.2909 3.9488 2.9649 
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Fig. 5. Relation between actual PCI and predicted PCI for all ANN models: a — model with 8 neurons; b — model with 9 neurons; 
c — model with 10 neurons; d — model with 11 neurons; e — model with 12 neurons; f— model with 13 neurons 


Looking at the results of the statistical analysis of the outputs of the developed models represented in Table 4 
and Figure 5 (a-f), we find that they have an acceptable ability to predict the condition of pavement based on the 
variables that were used in training the models. Also, notice that the best results were for two models with the 
structure (9-9-1) and (9-10—1), as they had the highest coefficient of multiple correlations (R°), which is 0.9493, 
0.9451, as well as the lowest values of MAE, RMSE, and MARE. 

Evaluation of pavement deterioration via ANN model: 

To recognize that the developed neural network models can help in predicting the rate of pavement deterioration, the 
idea arose to select several asphalt sectors and determine the effect of changing the operational life on the condition of 
each sector. To achieve this goal, each model defined the input values for each sector and fixed them, except for one 
variable, which was the operational life. Then we left the opportunity for the models to calculate the output values based 


on their previous experience. The sectors that were subjected to this experiment were numbers 504, 1227, 430, and 544, 
as shown in Figures. 6 (a, b, c, and d). 
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Fig. 6. Pavement sections deteriorations against service life: a — deterioration of section no. 504; b — deterioration of section 
no. 1227; c — deterioration of section no. 430; d — deterioration of section no. 544 


It is clear from the previous figures that the developed models have the ability to predict the condition of the asphalt 
as a result of the change in the service life. Where the value of the PCI reduces as the operational life increases, which is 
consistent with what happens in the field. 

This research may help those in charge of developing a road maintenance plan in order to limit the deterioration of 
the pavement condition by anticipating the pavement condition in advance and intervening promptly by choosing the 
appropriate treatment method in proportion to the capabilities available to the authority entrusted with this matter. 

Discussion and Conclusions. To achieve the objectives of this research, artificial neural network techniques were 
used to construct condition data deterioration models for pavement surfaces to assist decision-makers to choose the 
most efficient pavement maintenance solutions. Artificial neural networks (ANN) were utilized together with back- 
propagation algorithms to train the forecasting models for selecting appropriate options for treatment based on the value 
of the pavement condition index (PCI). A wide range of 1,614 cases of pavement condition description were provided 
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from previous investigations for various sectors of the Russian highway network to use as a database helping in creating 
the developed ANN models. Using MATLAB software, a large number of models with different numbers of neurons in 
the hidden layers were established, trained, and the results validated. The results of the models were studied, and the 
following conclusions were made: 

— ANN models showed their ability to represent the nonlinear relationship between a large number of variables and 
the pavement condition coefficient with high accuracy, as indicated by the high value of the multiple correlation 
coefficient (R’) for all stages of model learning (training, validation, and testing). 

— According to models’ validation, the variations between the actual PCI values and the predicted output results of 
the models were close to each other. 

— Evaluation of models generated with different numbers of neurons in the hidden layers showed that they all 
assumed close patterns as well as similar flows of surface deterioration rates. 

— The developed models save a lot of time and effort to determine the condition of the pavement and reduce the 
risks for the examiners compared to the visual examination method due to its dependence on the results of the dropped 
pregnancy test. 

— The obtained results suggested that decision-makers might use the created modeling to optimize maintenance or 
rehabilitation methods for the most impacted roads to lower the rate of degradation based on available resources. 
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Abstract 

Introduction. Determination of mechanical properties of layered structures of highways is an urgent task. This is due, 
firstly, to the need to control the quality of new sections during the construction of highways. Secondly, to assess the 
condition of existing roads with the accumulation of damage and defects is of interest. The formation of multiple 
defects (cracks) changes the averaged viscoelastic properties of the components of the structure, specifically, the surface 
asphalt-concrete layers. The article discusses the use of neural network technologies to improve the accuracy of the 
recovery of viscoelastic properties. This approach is based on experimental methods. As an example, we can give the 
definition of the dynamic deflection of a structure from a falling weight, FWD. 

Materials and Methods. The elastic modulus of a three-layer structure was determined on the basis of a neural network. 
To find out the solution accuracy, it was compared to the results of mathematical modeling and experimental data. 
Results. The experimental and calculated parameters of the elastic modulus of individual layers of the road structure 
turned out to be very close. The proposed approach to determining the mechanical properties of materials of road 
structures allowed us to apply the obtained results to examination of the condition of individual elements and the entire 
road structure. 

Discussion and Conclusions. The prospects of using artificial intelligence to determine the mechanical properties of 
layered structures was shown. Further improvement of methods and tools for analyzing the behavior of road structures 


under dynamic loading will expand existing approaches to assessing the condition of road structures. 


Keywords: neural network, non-destructive testing, modulus of elasticity, regression analysis, multilayer network, 


impact indentation, neural network technologies. 
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Introduction. Assessment of the state of road structures enables to predict the process of road operation, plan the 


timing and volume of repair and restoration activities, affect the efficiency and safety of road transport. Non-destructive 
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testing methods are actively used’. In a number of scientific papers [1, 2], the authors noted the reliability and cost- 
effectiveness of this approach. The non-destructive methods include the determination of the dynamic deflection of a 
structure from a falling weight (FWD). This approach involves the registration of the degree of deflection of the road 
structure surface at a certain distance from the impact point. At the same time, during the installation, the elastic 
modulus of the layers is calculated, the load and thickness of the layers of the structure are determined. 

One of the major scientific research objectives in this direction is to provide timely detection of defects and changes 
in the components of road structures. The researchers have identified the factor that influences the main characteristics 
of the dynamic stress-strain state of the “road structure — ground” system more than others. This is a change in the 
dynamic elastic modulus of the elements of the pavement, the base, and the ground [3-5]. 

An alternative tool for monitoring the state of road surfaces can be a neural network trained according to the results 
of experiments [6, 7]. It is worth noting that neural networks are extensively used in the construction, including road 
construction. Neural networks are an effective tool for solving some classification problems, including construction. 
They help to perform regression analysis as well. The traditional approach to the construction of nonlinear models has a 
number of limitations. They are connected, firstly, with complex, hidden nonlinear relationships of the source data. 
Secondly, over time, it is possible to modify the structure of the links, which necessitates the transformation of the 
nonlinear model structure. The use of neural network technologies makes it possible to eliminate the disadvantages of 
the traditional approach. A neural network is able to independently determine functional dependences, approximate 
nonlinear functions with a given accuracy, and rebuild its architecture in the learning process. The above verifies the 
timeliness of the presented work. 

Neural network technology opens up new opportunities for solving problems of design, optimization and diagnostics 
of building structures [8—11]. O. M. Maksimova writes that the possibility of additional training of the neural network 
provides refining the construction structure, improving its management [12]. In [13], it is noted that artificial 
intelligence can reduce the costs of companies, thanks to the timely diagnosis of construction projects. In [14], a broad 
overview of the possibilities of using neural networks in construction is presented. The author notes the multivariant use 
of artificial intelligence to solve a variety of complex nonlinear problems of forecasting, evaluation, optimization. 

In [15], special attention is paid to the quality and state of the road network. Based on experimental data, two 
models are compared: multidimensional analysis and artificial neural network. They were developed to create a 
decision support system for assessing the operational condition of roads due to the hydrogeological situation. The 
author concludes that the neural network model is more reliable and efficient. 

The selection of the type of neural network is of particular note. The authors [16] conclude that the type of neural 
network does not significantly affect the solution to the problem of calculating the elasticity modulus of the pavement 
layers. At the same time, the most similar results between the model and the target values were provided by a 
feedforward neural network model. 

The study aims at solving the problem of reverse calculation of the elastic moduli for the dynamic deflection 
bowls from the impact of the FWD unit. The results of simulation of dynamic impact on the road structure, and 
those obtained by a neural network, are compared. 

Materials and Methods. The dynamic effect of the FWD unit was simulated with the help of the program for 
calculating wave fields in the structure. The model of axisymmetric deformation of the medium used in the calculation 
can be analytically described as a multilayer half-space D, consisting of N layers in a cylindrical coordinate system 
(R,6,Z). The axis of symmetry OZ is directed vertically upwards and is orthogonal to the boundary plane of the half- 
space and its layers (Fig. 1). 


' Eremin RA, Kulizhnikov AM, Pudova NG. Metodika kompleksnogo obsledovaniya dorozhnoi odezhdy avtomobil'nykh dorog s primeneniem 
metodov nerazrushayushchego kontrolya. In: Proc. 16th All-Russian Sci.-Pract. Conf. of Survey Contractors. Moscow: Geomarketing; 2021. P. 438- 
446. (In Russ.) 

* Qinwu Xu. Development of a Computational Method for Inverting Dynamic Moduli of Multilayer Systems with Applications to Flexible Pavements 
under FWD Tests: dis. ... doc. of philos. Texas, 2014. 223 p. 
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Fig. 1. Multilayer construction of pavement in longitudinal section 


The medium motions are described by the dynamic Lamé equations in displacements for a homogeneous linear 
elastic material within each layer: 

(A + 2u) grad div u(r,t) — protrotu(r,t) = pii(r,t), (1) 
where wu = (Up, Uz) — displacement vector, r — radius vector of the observation point, A, — Lamé elastic 
coefficients, p — material density. 

The layers of half-space D are rigidly fastened to each other. When passing through flat interface boundaries, the 
displacement and stress vectors must match (0,,T,,-). 


Each layer of the half-space has the following characteristics: thickness h;, density ;, elasticity modulus E; = 
BAj+2 pj 
ante 


0,50; 
Athi 


, Poison's ratio v; = 


Boundary conditions for stress vector s = (0,,T,,) at the upper boundary of the multilayer medium D, distributed 
over the spot radius 79 of the contact between the wheel and the road structure, can be written as: 
. ttt 
Pysin (=),t E [0, c] . 
0,t>T 


where Py — constant intensity of impact loading with duration T. 


0, (R, t) = zr = 0, (2) 


To solve the problem, the integral Hankel transform is used in the radial Fourier coordinate and in time. As a result, 


a representation of the wave field is obtained: 
u(r, t) = fe dwexp(—iwt) i Q(u, w) + S(u, w)Jo1(uR)udu. (3) 

Here, S(u, w) — the Fourier—Hankel transform from the external load vector s; Q(u, w) — Green’s matrix—function 
for a multilayer half-space obtained through the superposition principle as a sum of matrices for homogeneous half- 
space, as in [4]. 

To analyze the properties of the structure in solving the inverse problem, the results of experimental measurements 
of the bowl of dynamic deflections at specified observation points on the surface R = d,, k= 1, 2, ..., 10, are used. The 
average distance between them is 0.3 m. The measured characteristic at each point on the surface: 

max |u, (dx, t)| 
. . 


Figure 2 shows an example of the behavior of vertical deflection functions over time for a typical three-layer road 
structure. 
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Fig. 2. Time dependences of vertical deflections uz(d,, t): a — simulation data; b — experimental data 


The procedure for restoring the elastic modulus of structure £; is based on minimizing the quadratic error function: 
FE) = Dida (Ee(Ed — 8?) 
where rh were obtained from an experiment using the FWD unit. The Levenberg-Marquardt method was chosen as 
the minimization algorithm. The elements of the gradient and the Hesse matrix of function F were calculated through 
finite differences of the Ist order. At the same time, the output of elastic modules E; beyond the physical ranges of 
changes in these parameters was controlled. 
To solve the problem with the help of neural networks, a standard pavement design was considered (Table 1). 


Table 1 


Experiment pavement design 


Structural elements Thickness, cm Piste odie eae es 
upper bottom 
Asphalt concrete layer 10 4,200 1,200 
Crushed-stone layer 44 1,500 100 
Ground 00 300 20 


Table 1 shows the main structural elements of the road section under study. The model was constructed with a 
maximum elastic modulus (4,200 MPa), since the experimental measurement results were tied to an average air 
temperature of 20 °C. To determine the elastic modulus in the winter period, it is required to recalculate the indicators. 

A fully connected neural network was implemented in Python using Keras deep learning library. Basic network 
parameters were as follows: 

— 300, 30, 3 neurons in layers; 

— 3,000 eras of learning. 

The activation function in the hidden layers was ReLu, in the output layer — sigmoid. The network processed 
normalized data in the range from 0 to 1. After network calculations, denormalization into a numerical form was 
performed. The neural network was trained on a data set of 1,000 values obtained under modeling. On average, the error 
of neural network calculations on the test set was within 10%. 

The neural network inputs were bowls of dynamic deflections. Figure 3 shows them at different moduli of elasticity 
of the layers of the road structure. The network outputs are values of elastic moduli of structural elements of the 


pavement. 
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Fig. 3. Bowls of dynamic deflections with different elastic modules 


Experimental data for the study were obtained using FWD unit on a road section with an asphalt concrete coating thickness of 
10 cm, and a crushed stone layer thickness of 44 cm. The elastic moduli of the layers of road clothes were determined (Table 1). 
The layout of the sensors-geophones: 0 — 20 — 30 — 45 — 60 — 90 — 120 — 150 — 180 — 210 cm from the impact 
point. 


Research Results. The authors have elucidated various combinations of structural elements: 

— asphalt concrete layer + crushed stone layer; 

— asphalt concrete layer + ground; 

— crushed stone layer + ground; 

— asphalt concrete layer + crushed stone layer + ground. 

Let us study visualizations of the results of neural network calculations. Figure 4 shows the results obtained by the 
neural network during the test calculation for the 1st option. The accuracy of calculations exceeded 90 %. 
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Fig. 4. Result of neural network training in option 1 — “asphalt concrete layer + crushed stone layer” 


As for the test calculations for options 2 and 3, in general, the calculation results were close to the first one. Figure 5 
shows a good result of neural network training on all layers of the structure. The presented values were obtained on the 
test set by a random sample. 
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Fig. 5. Error of neural network calculations under option (1) 289 
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The graphical representation of the ratio of the expected and the obtained result of the elastic modulus in Figure 5 
indicates insignificant errors in the test calculation of the layers of the road structure for option 1. 
Special attention should be paid to the text calculation for option 4. 
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Fig. 6. Result of neural network training in option 4 — “asphalt concrete layer + crushed stone layer + ground” 


Figure 7 shows the statistical distribution of the error of neural network calculations. Note that the average error of 
most measurements is about 10 %. However, in some cases, there is maximum deviation close to 23 % (on the 
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Fig. 7. Error of neural network calculations under option 4 — “asphalt concrete 
layer + crushed stone layer + ground”: a — at E2 € [1,200; 4,200] MPa; b — at E, € [100; 1,500] MPa; c — at Ep € [20,300] MPa 


From the data in Figure 7, conclusions can be drawn about the accuracy of neural network calculations on the test 
set. In the first case (Fig. 7 a), we note that of the entire test set, 41 % (82 measurements) had an error in the range of 
5-10 %. Similar indicators of 52 measurements (26 %) — 2-5 %. In 30 cases (15 %), values of 0-2 % were recorded. 
36 measurements (18 %) yielded results of 10-20 %. 

In Figure 7 b, the maximum indicator was 70 values for errors of 2-5 %. According to Figure 7 c, the results of 
calculations of the elastic modules of the ground had the smallest error. Thus, there are 116 values within 0O—2 %. This is 
58 % of the test set. 44 and 16 values are within 10 % of the error (22 % and 8 %, respectively). 

Discussion and Conclusions. Hence, neural network technology has proved to be an effective tool for determining 
the mechanical properties of materials of road structures. When solving the inverse problem, a high degree of recovery 
of the elastic modulus values was noted. The average error of calculations did not exceed 10 %. 

The neural network training options according to the values of the elastic moduli of the layers were studied. The 
error recorded in these cases was mainly due to insufficient training of the network. To solve the problem, it is required 
to increase the amount of data for training. 

The results of the neural network technology application were consistent with the experimental data obtained after 


the operation of the FWD dynamic loading unit. 
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